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ABSTRACT 


The motion of charged particles in a model magnetosphere is studied using 
the three adiabatic invariants. The particle shell geometry is determined, and 
drift velocities, bounce periods and equatorial pitch angles are computed as a 
function of local time. The following conclusions were reached: 

1. Shell splitting in the outer magnetosphere becomes important beyond 5 
earth radii: dipole-type descriptions of the radiation belt become invalid. 

2. Equatorial pitch angles tend to align along field lines on the night side of 
the magnetosphere, and perpendicularly to the field, on the day side. 

3. There are regions in the magnetosphere, where only pseudo-trapped 
particles can mirror, i.e. particles which will leave the magnetosphere before 
completing a 180° drift. 

4. Longitudinal drift velocities depart considerably from the dipole values 
beyond 5 R ( , . and can be as much as 2-8 times greater on the night side than on 
the day side. Thus a given particle spends 2-3 times more time in the day side 
than in the night side. 

5. The action of a pitch angle scattering mechanism will lead to a radial 
diffusion of particles. The loss mechanism will be greatly enhanced by scatter- 
ing of mirror points into the pseudo-trapping regions. 

6. After recovery from a prototype magnetic storm, particles which were 
in the day side during the sudden commencement will have higher energies, their 
shells having moved radially inwards. Particles caught in the night side will have 
moved outwards, with their energies decreased. 

7. The repeated action of magnetic storms will result in a net inward dif- 
fusion of particles, with a net increase of their energy. 

iii 



ON THE ADIABATIC MOTION OF ENERGETIC PARTICLES 


IN A MODEL MAGNETOSPHERE 


I. INTRODUCTION 

Recent experimental results on trapped particle flux behavior in the outer 
magnetosphere indicate that physical processes governing particle diffusion and 
acceleration, are strongly influenced by the trapping field itself and by the time- 
changes of its configuration. The main evidence comes from the observed strong 
correlations between particle flux and energy spectra variations beyond 2-3 
earth radii, with geomagnetic perturbations such as the sudden commencement 
of a geomagnetic storm or the ring current during the main phase [ Frank , 1966; 
Mcllwain , 1965; Mcllwain , 1966]. It seems therefore useful to attempt a detailed 
theoretical description of the behavior of a flux of trapped particles using a reason- 
ably accurate magnetospheric field model, and simulating prototype time varia- 
tions of the field configuration. The first detailed studies of this type were done 
by Hones [1963] for auroral particles, and by Fairfield [1964] for energetic 
particles. 

There are several sources of the field in the geomagnetic cavity: the mag- 
netization of the earth’s interior, the currents flowing on the surface of the mag- 
netopause, the currents in the "neutral sheet" of the tail of the magnetosphere and, 
eventually, diamagnetic ring currents originating in trapped particle density 
gradients at 2-4 earth radii. At geocentric distances of less than, say, 4 earth 
radii, the internal geomagnetic field dominates; beyond 4R e , the currents in 
the magnetopause (and in the neutral sheet) perturb the dipole-type internal 
field, and introduce a strong noon-midnight asymmetry. Any model must take 
these sources into account. 
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Before adopting de facto a given field model, let us list our requirements. 
First of all, we are mainly interested in particles trapped onfield lines reaching out 
beyond, say, 5 earth radii, on the equatorial plane. This means that we can safely 
ignore all higher multipoles of the internal field, and replace it by a centered 
dipole. Second, we shall also ignore the effect of a ring current. Third, we shall 
consider the dipole axis perpendicular to the sun-earth line, which of course is a 
very substantial limitation. However, a "wobbling" dipole would make our calcu- 
lations immensely more complicated, without, however, adding much to the gen- 
eral results, at least within the scope of this paper. Finally, electric fields will 
be ignored; this means that we are restricting ourselves to particles of high 
enough energy to ensure that the gradient drift always prevails over the E x B 
drift. 

A model which satisfies these requirements and which has already predicted 
or explained experimental results with good quantitative agreement, is that given 
by [ Mead , 1964; Williams and Mead , 1965; Mead , 1965] . This model considers 
two sources, in addition to the internal dipole; currents in the magnetopause, 
and currents in the tail of the magnetosphere. Four adjustable parameters de- 
termine the field in Mead's model; (1) The distance R s from the center of the 
earth to the magnetopause, in the solar direction; (2) and (3), the distances R m . n , 

R from the center of the earth to the close and far limit of the neutral sheet 

max 

in the anti- solar direction, respectively, and (4) the field intensity By near the 
neutral sheet. Most of the typical variations of these parameters. See Section IV 
for choice of parameters actually used. 

Fig. 1 shows field lines of this model in the noon-midnight meridians, cor- 
responding to the parameters which we shall adopt as describing the quiet-time 
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Figure 1-Closed field lines in the model magnetosphere, for parameter values given in the text (page 18). Field lines 
labeled with the geomagnetic latitude of their intersection with the earth’s surface. 




state of the magnetosphere. In this figure, all field lines start at equally spaced 
latitudes on the earth's surface. In order to save time and space, we shall call 
Mead's model of the magnetospheric field "the meadosphere." 

II. ADIABATIC INVARIANTS 

The motion of charged particles in a trapping field geometry can be described 
by means of three adiabatic constants of motion [Northrop, 1963]. Although the 
the invariants yield incomplete information about the actual position of a particle 
as a function of time, they do lead to the determination of the two dimensional 
manifolds or "particle shells," on which the guiding centers of the particles are 
confined. They further provide general information on the energy changes of the 
particles, although, again, no "microscopic" time-history would be available. 

The three adiabatic invariants are the following: (1) the magnetic moment 
M of the particle, generated by its cyclotron motion around a field line; (2) the 
second invariant J , associated with its bounce motion along a field line between 
mirror points; (3) the flux invariant, associated with its azimuthal or longitudinal 
drift motion. The definition of these quantities are: 


M = 


P[_ 

2m 0 B 


( 1 ) 



Pr, ds 

A • dx 


( 2 ) 

(3) 
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and p n are components of the momentum perpendicular and parallel to the 
magnetic field vector, respectively; B is the absolute magnetic field at the 
instantaneous position of the guiding center and m 0 the rest mass. In (2), the 
integration is extended along the field line for a complete bounce oscillation; 
ds is the element of arc of the field line. (3) represents the magnetic flux en- 
closed by the particle shell; X is the magnetic vector potential, and the integral 
is extended along any closed path passing around the particle shell and lying in 
it. These quantities are adiabatic constants, i.e. conserved only under certain 
conditions. Each of the invariants has an associated characteristic period of 
time and a characteristic length. These are respectively; (1) cyclotron period 
and gyroradius; (2) bounce period and arc length between conjugate mirror points; 
(3) azimuthal drift period and arc length of the equatorial ring of the shell. In 
most magnetic field configurations, these characteristic quantities differ by sev- 
eral orders of magnitude from one another. Adiabaticity requires that the field 
configuration should not change appreciably during a characteristic period. If 
this condition is violated, the corresponding value of the invariant will no longer 
be conserved. The other two may still remain unaffected. 

In place of (1) and (2) we shall introduce two other expressions which are 
much more convenient for numerical computations. According to (1), we intro- 
duce the mirror point field intensity B m : 


B 


m 


2m 0 M 


(4) 


p is the total momentum of the particle. Further, we introduce the geometric 
integral 
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(5) 


I = 



*L>ds = f 

B 2 P 


extended along the field line between one mirror point and its conjugate. B m and 
I uniquely determine a particle shell and will be used in what follows as the two 
identifying parameters for such a shell. Notice at once that for a time-independent 
magnetic field, in absence of electric fields, (4) and (5) are adiabatic invariants, 
too. The advantage of (4) and (5) is that they only depend on the field geometry. 

For time dependent fields, care has to be taken regarding relation (5). If p 
changes appreciably during one bounce, we have: 


= 4r f P(s) 1 /l~- 

2p 2p / V 


B(s) 

B 


ds / I 


Only for slowly varing fields, (5) holds at all times. 


(5') 


Introducing the usual relativistic factor y - m /m 0 , we can combine (4) and 
(5) with (1) and (2) to obtain 


(y 2 - 1) I 2 = Kj = const . 

I 2 B = K , = const . 

m 2 

To these, we must add (3), which we write in the form: 


( 6 ) 

(7) 


<j> = $ (B , I) = K 3 = const . ( 8 ) 

These are the adiabatic constants which will be used henceforth. They uniquely 
determine y, I and B m for any static field configuration. For time dependent 
fields, Kj and K 2 in general are not constant during the .'nterval of change (see 
(5*)). However, if this interval is transient, i.e. if it is preceded and followed by 
time independent states of the field, the constancy of K and K 2 does hold when 
(6) and (7) are evaluated for the initial and the final states. 
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For a field constant in time, (8) is no longer needed, and y, I and EFare 
conserved individually. A special case is that of particles mirroring on or close 
to the equator (I = o, B m = B e )• In that case we replace (6) by the following, 
derived from (4): 


-* B = const. 

y 2 - 1 e 


(6a) 


This holds all the time, even in non-static fields. In this special case, (8) is a 
function of B e only. 

We have to add Liouville's Theorem, to relations (6), (7) and (8), in order to 
complete our description of trapped particle dynamics. We introduce the direc- 
tional, differential flux of particles populating a given I , B m shell: 

(9) 

j = j (7» I. B m ) Part icles/st e radian x energy x sec 

Assume particles distributed on shell in such a way that there is a steady state. 
For each point of the shell, there is a unique cone along whose elements the given 
group of particles is streaming. This direction is the particle's pitch angle, which 
at the equatorial points of the shell is given by 



3 e = B e (I, B m ; azimuth) is the minimum or equatorial B value of a particular 
field line of the shell. We must point out at once that even for time-independent 
fields, a g is not an adiabatic invariant like the mirror point field B m ; in general 
it will depend on azimuth (longitude or local time) through B e . 

The particle density in phase space, f , associated with the flux (9) is given by 
f = j/p 2 . Liouville's Theorem states that this density in phase space remains 
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constant along the dynamical path of a particle. This means that, as long as all 
particles always remain on a common shell*, even if the shell itself changes with 
time, the following conservation theorem holds: 


iSll ll B m) = K . = const . 
y 2 - 1 4 


( 11 ) 


(6), (7), (8), and (11) are the basic expressions to be used for the study of the 
evolution in space and time of trapped particles in the outer magnetosphere. 

The first three lead to the determination of the actual shell of a given group of 
particles, and their actual energy; equation (11) gives the actual value of the 
directional differential flux of this group of particles. 

As was shown by [ Anderson , Crane , Francis , Newkirk and Walt , 1964] the 
average equatorial azimuthal drift velocity u of the particles populating and I, B m 
shell can be obtained as a direct consequence of Liouville’s Theorem. The 
value is, with our notation: 


m Q c 2 y 2 - 1 VI 
~eB y ST 

e / b 


( 12 ) 


VI is the limit of SI/Sy where 8y is the equatorial distance between two neigh- 
boring shells, each one characterized by I and I +8 1, respectively, and by the 
same B m -value. S b is the half-bounce path, i.e. the rectified path of the particle 
between one mirror point and its conjugate: 


S 


b 




*This condition will be fulfilled if all invariants (6) — (8) are conserved throughout the evolution 
of the system. 
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S b is related to the bounce period by r b = 2 S b /v where v = particle velocity. 
We can obtain an expression for S b by taking the derivative of (5) with respect to 
the mirror point field intensity B m along a given field line. It can be shown by 
simple algebra that 


S 


b 


I - 2 B 

m 


BI 

3B 


(13) 


The derivative has to be taken along the given field line. This expression has 
general validity for any trapping field geometry, and is very useful for computa- 
tional purposes, for it only requires the calculation of I on two neighboring 
points of a field line. 


ID. QUALITATIVE DISCUSSION 

Before getting into the discussion of numerical computations for particle 
shell geometry and time variations of the meadosphere, it is useful to present 
a qualitative analysis of the consequences of the previous section. First of all, 
let us envisage a trapping magnetic field constant in time. In this case, (4) and 
(5) are conserved. We can assign to each point in space a pair of values I, B m 
such that a particle mirroring there has the value I for the integral (5), B m 
being simply the field intensity at that point. In this way, I and B m become uni- 
form functions of space. As the particle drifts to other field lines, it must keep 
these values constant, i.e., it will cover a shell of field lines which pass through 



Figure 2-Geometric definition 
of o particle shell 
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the intersections of two given constant-I and constant-B m surfaces (Fig. 2). No- 
tice carefully that the surface I = const, is not the particle shell. 

Let us consider the geomagnetic field. Take a particle which starts at a 
given longitude <p , circling around a given field line and mirroring at a value B m . 
The integral (5) computed along the field line between the two mirror points, has 
a value I. This means that when drifting through any other longitude, say 180° 
away, this particle will be bouncing along a field line which passes through the 
intersection of the corresponding I = const, and B m = const, surfaces. Now take 
a particle which starts on the same initial field line, but which mirrors at a lower 
value B' m <B m (Fig. 3). Its integral (5) will also be smaller, I' < I. After a 180° 


I’n CONST 



AT 4 > 



SPLIT SHELLS 
AT 0 + 180° 


Figure 3-Shell splitting in 
asymmetric fields. 


longitudinal drift, this second particle will be traveling along a field line which 
passes through the intersection of the surface I' = const, and B' m = const . 

Only in case of perfect azimuthal symmetry (as in the pure dipole), will these 
surfaces intersect exactly on the same line as that of the first particle. In the 
general case, particles starting on the same field line at a given longitude will 
populate different shells, according ta their initial mirror point fields, or, what 
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is equivalent, according to their initial equatorial pitch angles (10) (of course, all 
these different shells would be tangent to each other at the initial field line). 

For the case of the real geomagnetic field in absence of external perturba- 
tions, i.e. within about 3-4 earth radii, it can be shown that the distance between 
split shells is only very small, a fraction of 1% of the distance of the equatorial 
point of a field line to the center of the earth. In other words, with a very good 
approximation, one can say that all particles initially on the same field line, will 
mirror on a common field line at any other longitude. This has an important 
consequence: it enables a two-dimensional description of the three-dimensional 
radiation belts, at least up to distances of about 4R e . Indeed, if particles do popu- 
late the same shell irrespective of their initial mirror points, the omnidirectional 
flux of these particles will be the same on all points of the shell having the same 
B value, irrespective of the longitude (i.e. local time) (provided of course, no ap- 
preciable injections or losses occur during the drift). In order to describe 
omnidirectional particle fluxes in the inner magnetosphere, we therefore need 
only two "space" parameters: the value of the magnetic field intensity at the point 
of measurement and a parameter which characterizes the (unique) shell which 
goes through that point. This latter is Mcllwain's L-parameter [ Mcllwain , 1961]. 

L is a particular relation between I and B m which remains constant (within < 1%) 
on a given field line, and, therfore, on the whole shell generated by particles start- 
ing on that field line. Numerically, L gives the average distance of the equatorial 
points of a shell to the magnetic center. 

But what happens in the outer magnetosphere, where the azimuthal symmetry 
is brutally removed? Particles starting on the same field line, say in the noon 
meridional plane, will now populate different shells, depending on their initial 
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mirror points or equatorial pitch angles. For instance, they will cross the mid- 
night meridian on different lines. 

Let us start with a particle mirroring at or near the equator, on a line in 
the noon meridian, close to the boundary. For this particle, I L 0; according to 
(6a) it will drift around the earth on the equator following a constant- B path. 

This constant- B path comes considerably closer to the earth at the night side, 
because the field is weaker there (less compression), and we must go to lower 
altitudes in order to find a given B value. On the other hand, a particle which 
starts on the same field line on the noon meridian, but which is mirroring at 
high latitudes, will have a high I value. Under these circumstances. Mead has shown 
[ Mead , 1965] that the value of I is not much different from the arc length of the field 
line between mirror points. On the midnight meridian, the particle will there- 
fore be found on a line which has nearly the same length than the initial one, i.e. 
stretching out to roughly the same equatorial distance. In summary, all particles 
initially on the same noon-line, will cross the midnight plane on line portions 
sketched in Fig. 4. Furthermore, it is easy to realize that particles mirroring 


N 



Figure 4-Qualitative picture of shell 
splitting in a model magnetosphere for 
particles starting on a common field 
line in the noon meridian. 
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inside that area (BB' ), will cross the noon meridian outside (AA') of the initial 
line. If this noon-line is very close to the boundary, no stably trapped particle 
could be found mirroring inside the hatched area in the midnight meridional 
plane. Any particle doing this would not be able to complete a drift around the 
earth: it would abandon the magnetosphere before reaching the noon meridian. 
We shall call this a "pseudo-trapped" particle 11 (only transiently trapped). No- 
tice finally that a sharp trapping boundary in the noon side would not result in a 
sharp boundary in the back side. 

On the other hand, for a given field line in the midnight meridian, all par- 
ticles mirroring anywhere on this line, will cross the noon meridian in an area 
like the one shown in Fig. 5. All particles mirroring outside that area (BB' ) 
will cross the midnight-meridian outside (AA' ) of the given line. If now there 

N 

Figure 5 — Qualitative picture of 
shell splitting for particles start- 
ing on a common line in the mid- 
night meridian. 



is an "obstacle" behind that line (like for instance the neutral sheet), no stably 
trapped particle could be found outside the hatched area in the noon meridian. 
Any particle injected there, would be lost into the "obstacle" before reaching 
the midnight meridian: in this high latitude noon region, only pseudo-trapped 
particles could exist. 
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We can also make some qualitative predictions regarding equatorial pitch 
angles and longitudinal drift velocities of particles trapped in the meadosphere. 
Take again a group of particles starting on a common noon-side line. It is 
easy to see from Fig. 4, that the equatorial B-value of a given shell must 
decrease towards the midnight meridian except for particles mirroring close 
to the equator. This means that, according to (10) , the equatorial pitch angle 
must also decrease: particles align along field lines, in the midnight meridian. 
On the other hand, inspection of Fig. 5 reveals that the opposite is true for 
particles starting on a common field line in the midnight meridian: they 

align perpendicularly to the field line on the day side of the meadosphere. 

Regarding the drift velocity (12), we can qualitatively say that the field near 
the day-side boundary is much more homogeneous than on the night side, or than 
in the case of a pure dipole (Fig. 1). Therefore, VI will be relatively smaller, 
and S b greater. Both facts indicate that the particles will drift slowest on the 
noon meridian, spending therefore a greater fractional time in the day side of 
the meadosphere. This may have very important consequences for outer belt 
dynamics. 

Turning now to a qualitative discussion of the effects of a time-dependent 
magnetic field, we first have to point out that slow changes, which conserve 
all three invariants (6), (7) and (8), are reversible. This means that whenever 
the field is back to its initial configuration, all particles will be back on their 
initial shells with no net change in energy or directional flux. During the 
change, itself, a given shell will be deformed, its particles being accelerated 
or decelerated and their directional flux changing accordingly. In order to have 
a net change in shell, energy and flux after a complete cycle, it is necessary 
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to have a non-adiabatic change in the field, which violates, say, the third invari- 
ant. In this case, the time history of a particle depends on where in longitudinal 
position on the shell the particles was surprised by the non-adiabatic field 
change. It is easy to realize that due to this, the end-effect will always be a 
diffusion, even if the cause itself (the non-adiabatic change in field configura- 
tion) is not at all a stochastic process. Non-adiabatic compressions and dis- 
tentions of the magnetosphere very likely are the main cause for radial diffusion 
and acceleration of trapped protons [ Nakada , Dungey and Hess , 1965]. 

There is, however, another possible radial diffusion process, associated 
with shell splitting in the meadosphere. Suppose the action of an elastic pitch 
angle scattering process, such as interactions with electromagnetic or hydro- 
magnetic waves. This is a very short-time process which violates all three 
invariants. However, the fact that the particle remains on the same field line 
during the event (within one cyclotron radius) predetermines the new shell on 
which the particle will drift after the interaction. Inspection of Figs. 4 and 5 
reveals that for instance, a scattering process which, occuring on the noon 
side, increases the pitch angle (lowers B m ), will bring the particle to a shell 
which on the average gets closer to the earth. The same scattering process, 
occuring on the night side, would situate the particle on a shell extending 
further out. Any type of pitch angle diffusion process will therefore be accom- 
panied by a radial diffusion, which will be inwards or outwards according to 
where in longitude the particles are more likely caught by the individual proc- 
esses. If the original pitch angle scattering process is elastic, there would be 
no change in energy in this type of radial diffusion. Such a diffusion mechanism 
was experimentally observed in the Laboratory, and studied theoretically, by 
[ Gibson, Jordan and Lauer, 1963]. 
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IV. NUMERICAL RESULTS 


A computer code was set up to determine particle shells in the meadosphere, 
in order to study shell splitting and longitude dependence of drift velocities and 
equatorial pitch angles for the static case, and in order to analyze the evolution 
of a system of particle shells in a time- dependent case. The computer program 
consists of four main parts, which perform more or less independent operations. 

I. Field line geometry. This part furnishes complete geometric informa- 
tion about a given field line, computes I and B m valuesfora set of particles 
mirroring on this field line with prefixed equatorial pitch angles [re- 
lations (10), (4) and (5)J and determines their drift velocities and 
bounce periods [relations (12) and (13) j . 

II. Shell geometry. Given the particles defined in part (1), this section of 
the program finds the points of prefixed I, B m values at other longi- 
tudes, and traces the corresponding field lines, computing for each 
case, drift velocity and bounce period. In this way, the complete shell 
for each particle is generated. 

HI. Third invariant. This part computes the magnetic flux (3) enclosed 
1 

by the shell generated by a given particle. It is a combination of (1) 
and (2), and a subroutine which computes the integral (3) along the 
equatorial ring of the shell. 

IV. Non-adiabatic compression. Given a particle on a field line, this 
part changes the field configuration simulating a sudden commence- 
ment, finds the new field line (going through the same intersection 
with the earth’s surface), and locates the particle’s mirror points 
supposing conservation of the first two invariants. 
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V. Shell deformation during adiabatic time variations. This part finds the 


new I and B m values, as well as the energy, of a particle after a slow 
adiabatic change in the trapping magnetic field. Conservation theorems 
(6), (7) and (8) are used in this computation. 

Part (1) is mainly based on Mcllwain's INVAR code, conveniently imple- 
mented for our purposes. Part (2) contains a key subroutine called SEARCH, 
which at a given longitude finds the point of prescribed I and B m values (within 
prefixed tolerances), by an iteration method. This subroutine is quite fast (a 
fraction of a second on an IBM 7094, for a relative tolerance of 10' 4 in I and B). 
Part (5) contains subroutine LOOK, which starts with an approximate value of I , 
finds B m through (7) for prescribed values of K t and K 2 (determined by the initial 
state) , and computes the flux (8). The value of I is then corrected and the proce- 
dure iterated until the prefixed value K 3 for $ is approached within the wanted 
tolerance; y is then computed through (6). Remember that each step, i.e. each 
evaluation of flux requires the complete determination of a particle shell. Even 
so, this program is quite fast (10-20 seconds, to find a prefixed K x , K 2 , K 3 shell) 
with a relative tolerance for 4> of 10* 3 ). 

This program was applied using Mead's model of the magnetosphere. The 
numerical values of the four intervening parameters (Section I) for the quiet 
meadosphere were taken from [Ness and Williams, 1966]: R = 10R ,R . =8R , 

R max = 200 R e andB T = 15 gammas. Typical closed field lines are shown in 
Fig. 1. The kink of the field lines reaching out beyond 8 R e in the night side is 
caused by the assumption in Mead's model of a two dimensional neutral sheet 
(recent measurements, however, suggest a finite thickness of several earth 
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EARTH RADII (R e ) 


Figure 6-Relationship between geomagnetic latitude of the intersection of a field 
line with the earth’s surface, and the radial distance to its equatorial point, for 
the noon, midnight and dawn/dusk meridians. 

radii [ Bame , Esbridge , Felthauser , Olson and Strong, 1966] . Figure 6 shows the 
relationship between geomagnetic latitude of the intersection of a field line with 
the earth’s surface, and the radial distance to its equatorial point, for the noon, 
midnight and dawn/dusk meridians. 

Parts I and II of the code were applied to this quiet time field configura- 
tion, to obtain magnetic shells for particles initially mirroring on a common 
field line, and having equatorial pitch angles with cosines 0.2, 0.4, 0.6, 0.8 and 
nearly 1 (mirroring close to the earth’s surface). Fig. 7 shows how particles, 
starting on a common line in the noon meridian, do indeed drift on different 
shells, which intersect the midnight meridian along the field lines shown in the 
figure. The dots represent particles’ mirror points. Curves giving the position 
of mirror points for constant equatorial pitch angles are traced for comparison 
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Figure 7— Computed shell splitting for particles starting on common field lines in the noon meridian. Dots represent 
particles' mirror points. Curves giving the position of mirror points for constant equatorial pitch angle a are shown. 








(in a dipole field, they are constant latitude lines). Notice the change (decrease) 
in equatorial pitch angle for the same particle, when it drifts from noon to mid- 
night. Fig. 7 clearly confirms our qualitative predictions given in Section II: 
shell splitting becomes considerable beyond 5 earth radii, and completely in- 
validates the use of "L- values" or any other dipole-type description of the 
outer radiation belt. Fig. 8 depicts the same features for particles starting on a 
common field line in the midnight meridian. In this case, again, the pitch angle 
changes considerably when the particle drifts to the opposite meridian (increas- 
ing at noon). 

Notice from Figs. 7 and 8 that as equatorial pitch angles increase, shell 
splitting is directed radially inwards for particles starting on the same field 
line at noon, and radially outwards for particles starting on a common field line 
at midnight. Furthermore, shell splitting is maximum for particles mirroring 
close to the equator: for a given change in pitch angle (in degrees), the radial 
displacement of the shell will be greater for equatorial particles. 

When the initial field line has its equatorial point beyond about 8 a frac- 
tion of the particles mirroring on it can only be pseudo-trapped, being lost before 
completing a 180° drift. In particular, the computations reveal that particles 
mirroring at low latitudes in the back side, abandon the meadosphere through the 
boundary 30-40 degrees before reaching the noon meridian. On the other hand, 
particles mirroring at high latitudes on the day side, run into the tail (open field 
lines) 10-20 degrees before reaching the midnight meridian. Fig. 9 shows 
computed limits between stable trapping and pseudo-trapping regions in the 
meadosphere, on the noon-midnight plane. At other local times, both regions 
approach more closely the meadospheric boundary; going from noon to midnight, 
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Figure 8-Computed shell splitting for particles starting on common field lines in the midnight meridian. Dots represent 
particles' mirror points. Curves giving the position of mirror points for constant equatorial pitch angle a are shown. 




60 ° 


90 ° 


60 ° 



Figure 9— Location of the “pseudo-trapping” regions in the magnetosphere. Particles 
mirroring inside those regions are unable to complete a 180 drift around the earth. 
Those injected into the left side will be lost into the tail; those injected into the right 
portion will abandon the magnetosphere through the boundary, on the day side. 
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one disappears at the expense of the growth of the other. These results indicate 
the existence of a quite considerable loss cone in the day side of the meado- 
sphere. 

Figs. 10a and 10b summarize the information about the shell splitting effect. 
In both figures, the relation between noon and midnight radial distances to the 
equatorial points of a particle shell is given. In Fig. 10a, particles start on a 
common field line at noon, reaching out to R noon earth radii; in Fig. 10b, par- 
ticles start on a line at midnight, reaching out to R midn . . Curves are labeled 
with the cosines of the initial pitch angles. Figs. 11a and lib show how these 
pitch angles change when the particles drift to the opposite meridian. Notice 
again the marked tendency of particles to align along field lines on the night side, 
and to "squeeze" transverse to the field on the day side. 

The numerical calculations also confirm our predictions for local time 
dependence of the equatorial drift velocity (Section II). The geometrical factors 
appearing in (12) were computed, and represented in Figs. 12a and 12b as a 
function of the equatorial distance of the corresponding field line, for different 
pitch angles, and for noon and midnight, respectively. For a better understand- 
ing, angular drift factors are shown. For radial distances < 3R e , we observe a 
dipole- like dependence. Beyond 3 R e , there is a considerable departure. Drift 
velocities on the night side are indeed appreciably higher than on the day side. 
The peculiar inversion of the pitch angle dependence, occuring on the day side 
(Fig. 12a), is due to a shift in relative importance of curvature drift versus 
gradient drift (particles mirroring at high latitude experience more of a dipole- 
type field, and drift faster) . 

In Figs. 13a and 13b we have represented the percentage change of the linear 

drift velocity of a given particle, when it drifts to the opposite meridian. A 
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Figure 1 la-Particles starting in the noon meridian at a field line reaching out to R and 
initially having equatorial pitch angles of cosines 0.2, 0.4, 0.6, 0.8 and 1.0 respectively, will 
appear on the midnight meridian with cosines of pitch angles given by the curves. Notice the 
effect of alignment along the field lines (occuring at midnight), for radial distances > 6. 



Figure 1 Ib-Particles starting in the midnight meridian at a field line reaching out to R noon 
and initially having equatorial pitch angles of cosines 0.2, 0.4, 0.6, 0.8 and 1.0 respectively, 
will appear on the midnight meridian with cosines of pitch angles given by the curves. Notice 
the effect of alignment along the field lines (occurring at midnight), for radial distances ^ 5. 








close inspection of the local time dependence of the drift velocity (not shown 
here) leads to the important conclusion that a given particle trapped in the outer 
meadosphere (/,6R e ) spends up to 2/3 -3/4 of its total lifetime in the day side. 

In other words, there is always a higher probability to find a particle in the day 
side than in the night side. It can be shown that, as a consequence, particle volume 
densities can be about two times greater on the noon meridian, than at midnight, for 
a given class of particles. This represents an additional important asymmetry for 
trapped particle fluxes in the outer meadosphere. 

As discussed in Section II, any pitch angle scattering mechanism will lead 
to radial diffusion, due to shell splitting. According to the preceding results, 
this radial diffusion will be most effective for equatorial particles. This mechan- 
ism would tend to mix and blurr energy spectra of particles at different radial 
distances. Furthermore, the existence of large pseudo- trapping regions, es- 
pecially on the day side, implies an efficient particle sink for any pitch angle 
scattering mechanism (enhanced loss cone). On the other hand, the reverse could 
also be true: particles which happened to enter the pseudo-trapping regions from 
outside, could be scattered into stably trapped orbits by any pitch angle scatter- 
ing mechanism. 

We now turn to the numerical results for a time-dependent magnetic field 
configuration. The purpose is to study the trapped particle behaviour during a 
simulated magnetic storm. Following steps were adopted as a "model" storm: 

1. sudden, non- adiabatic compression, simulated by a decrease AR s of the 
parameter R s (inward displacement of the magnetopause). 

2. optional sudden increase AB T of the tail field. 

3. gradual, adiabatic recovery to the initial field configuration. 
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For (1) and (2) it is assumed that particles, while violating the third invariant (8), 
stick to their original field line driven by the dominating ExB drift, still conserv- 
ing the first two invariants (this "original" field line is supposed to be rigidly 
rooted in the ionosphere, during the sudden compression). Part (4) of the 
computer code is applied for this calculation. The gradual recovery in (3) is as- 
sumed to be flux-conserving; part (5) of the computer code was used here. 

Some of the results are summarized in Fig. 14, in which the percentage change 
of kinetic energy of a particle is represented as a function of equatorial distance 
of the initial field line, for different pitch angles, and for particles caught by the 
compression at noon (upper curves), and at midnight (lower curves). This energy 
variation is independent of the initial energy, and increases with the amount of 
compression AR S (taken = 2 R e for the curves in Fig. 14). An inward motion of 
the boundary of only 1 R g would yield energy changes roughly 0.45 times those 
shown in the figure. The curves given in Fig. 14 were derived for a sudden com- 
mencement without a sudden increase of the tail field. If one adds a typical 
increase AB T of 15 gammas, acceleration, deceleration and radial displacements 
become greater by a factor of about 2, the effect being considerably enhanced for 
particles which were in the night side during the sudden commencement. 

One clearly sees in Fig. 14 that the final effect of a storm depends on where 
in local time the particle was caught during the non-adiabatic phase: particles 
which happened to be in the day side of the meadosphere will have a higher en- 
ergy when the field recovers to the initial state; particles surprised in the night 
side will be decelerated. Furthermore, the first group of particles will have 
moved radially inwards, whereas the second group will be found on shells "in- 
flated" outwards (Fig. 15). It should be pointed out that during the first phase of 
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PERCENT KINETIC ENERGY CHANGE (AE) 


PARTICLES 
STARTING 
AT NOON 


PARTICLES 
STARTING AT 
MIDNIGHT 


cosa, 


EARTH RADII (R e ) 

Figure 1 4— Percentage kinetic energy change after a prototype storm, for particles caught by the 
sudden commencement in the noon meridian, and at midnight, respectively, as a function of the 
radial distance to the equatorial point of the initial field line, and for different initial pitch 
angles. The inward displacement of the boundary was taken as A R = 2 earth radii. 




Figure 15— Relation between initial and final position of the equatorial point of a 
shell at a given meridian (noon and midnight), after a prototype storm, for par- 
ticles caught by the sudden commencement at noon, and at midnight, respectively 

(AR s .= 2). 
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compression, both groups of particles attain higher energies; it is during the 
adiabatic recovery, that this peculiar asymmetry arises. The fact that a given 
particle always spends more time in the day side, leads to the conclusion that a 
magnetic storm should always have a net effect of inward diffusion and acceler- 
ation of trapped particles, after the field has recovered to the same initial con- 
figuration. We finally must mention that during a magnetic storm a considerable 
fraction of particles can be driven into pseudo-trapping regions of the meado- 
sphere, and therefore be lost through the boundary, or into the tail. This effect is 
particularly important for storms with increases of the tail field, which mainly 
leads to losses through the boundary on the day side. Likewise, particles which 
for some reason happened to be injected into pseudo-trapping regions during re- 
covery, can become stably trapped, under favorable circumstances of injection. 

The repeated action of magnetic storms will therefore cause acceleration 
and radial diffusion towards lower altitudes, of particles trapped in the outer 
meadosphere. This mechanism leads to energy spectra which depend on radial 
distance, hardening towards lower R values. This is indeed observed for pro- 
tons [ Davis and Williamson, 1963; Vernov , Vakulov , Kuznetsov , Logatchev , 
Nikolaev , Sosnovets and Stolpovsky , 1966], and was explained theoretically by 
[ Nakada , Dungey and Hess, 1965]. Electrons, however, should in addition be 
subject to pitch angle scattering by electromagnetic waves; as discussed above, 
the radial diffusion which accompanies pitch angle scattering when shell splitting 
is considerable, should greatly blurr the radial dependence of energy spectra 
for electrons. 

IV. CONCLUSIONS 

1. Shell splitting in the outer magnetosphere becomes important beyond 5 
earth radii; dipole-type descriptions of the radiation belt become invalid. 
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2. Equatorial pitch angles tend to align along field lines on the night side 
of the magnetosphere, and perpendicularly to the field, on the day side. 

3. There are regions in the magnetosphere, where only pseudo-trapped 
particles can mirror, i.e. particles which will leave the magnetosphere before 
completing a 180° drift. 

4. Longitudinal drift velocities depart considerably from the dipole values 
beyond 5 R e , and can be as much as 2-3 times greater on the night side than on 
the day side. Thus a given particle spends 2-3 times more time in the day side 
than in the night side. 

5. The action of a pitch angle scattering mechanism will lead to a radial 
diffusion of particles. The loss mechanism will be greatly enhanced by scat- 
tering of mirror points into the pseudo-trapping regions. 

6. After recovery from a prototype magnetic storm, particles which were 
in the day side during the sudden commencement will have higher energies, their 
shells having moved radially inwards. Particles caught in the night side will 
have moved outwards, with their energies decreased. 

7. The repeated action of magnetic storms will result in a net inward dif- 
fusion of particles, with a net increase of their energy. 
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APPENDIX I 

Listing of the program briefly described in the text. For more 
details, see the comment cards and the list of output which follows 
the listing. The deck is available on request. 
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LIST OF OUTPUT 


Initial field line: 

VN 1: Radial distance in earth radii to field line points 

OLAT: latitude of field line points 

OLONG: longitude of field line points 

B: field intensity in gauss, at field line points 

SPITCH: cosines of equatorial pitch angles for a group of particles on the field 
line 

BDIM: mirror point field intensities for these particles 

FIDIM: corresponding values of second invariant (in earth radii) 

SALT , SLAT, SLONG: coordinates of mirror points 

SPATH: half-bounce path (expression 13 in the text) corresponding to these 

particles (in earth radii) 

SDRIFT: geometric factors giving drift velocities of these particles (for con- 
version into cm/sec, see comment card) 

EALT: radial distance to the equatorial point of the field line 

EQB: equatorial B-value 

ENERGY: kinetic energy of the particles, in units of rest energy 
FLX: magnetic flux subtended by the shell generated by the particles (in 

gauss (earth radii) 2 ) 

New field line (at a different longitude (splitting) , or at the same longitude after 
the simulation of a storm). 

VN1, OLAT, OLONG, B: same as above 

RPITCH, RALT, RLAT, RLONG, RPATH, RDRIFT; same as SPITCH, SALT, etc. 
EQALT: radial distance to the equatorial points 
EQBB: equatorial B-values 
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SPLIT 

COMMON B I 200 ) , VNl C200 ) v VN2( 200 ) , VN3 ( 200 ) , ARC 1 200 ) , VNE Aft ( 3 ) ,VB0(3) » A 2 

IVSAVE(3),B0,BNEAR, JUP,MMM A„ 3 

COMMON FRONT, SHEET I , SHEET2, FSHEET A A 

DIMENSION V ( 3 ) * BDIM(9>, FIDIM(9), EOALT ( 9 ) , EQLONG ( 9 1 A 5 

DIMENSION S P I TCH ( 9 ) f RPITCH(9) A 6 

DIMENSION SALT(9), SL AT ( 9 ) , SLONGI 9 ) , R ALT t 9 ) # RLAT19), RLCNG (9 ) A 7 

DIMENSION SDR I FT ( 9 ) , RDR I FT t 9 ) , SPAtH(9), RP ATH { 9 ) , EQBB (9) , RLOST A 8 

IT 19) A 9 

DIMENSION ENERGY ( 9 ) , FL X ( 9 ) » QV2I9), QV3I9) A 10 

DIMENSION BEG ( 20C ) r BEND1200), EC0(200), BL0GI200) A 11 

ERR =0.0001 A 12 

ERR B = 0 • 004 A 13 

ERR 1=9 • 008 A 14 

E R R F = 3 . 0 0 1 A 15 

A 16 

FRONT = DISTANCE TO MAGNETOPAUSE AT STAGNATION POINT IN EARTH RADI A 17 

SHEHT1 = MINIMUM DISTANCE TO NEUTRAL SHEET IN EARTH PADII A 18 

SHEET 2 = MAXIMUM RADIAL EXTENSION OF NEUTRAL SHEET IN MIDNIGHT A 19 

MERIDIAN IN EARTH RADII A 20 

FSHEET = FIELD STRENGTH PARALLEL TO NEUTRAL SHEET IN GAMMA A 21 

A 22 

READ (3,39) FRONT , SHEET1 , SHE ET2 , F SHEE T A 23 

A 24 

KM AX = NUMBER CF MIRROR POINTS AND PITCH ANGLES WANTED (.LE.9) A 25 

DELTA = STEP SIZE IN LONGITUDE (DEGREES) , A 26 

START = LONGITUDE INTERVAL TO BE SKIPPED BEFORE STARTING STEPS A 27 

TERM IS LONGITUDE WHERE SHELL TRACING SHOULD BE TERMINATED A 28 

A 29 

READ (5,40) KMAX, DELTA, START,TERM A 30 

A 31 

I ST CRM = 1 MEANS THAT AFTER FIRST SPLIT ANALYSIS, A GEOMAGNETIC A 32 

STCRM SHOULD BE TURNED ON. WIND IS THE INWARD DISPLACEMENT A 33 

OF THE FRONT SICE OF THE MAGNETOSPHERE. STORM IS THE I NWARD A 34 

DISPLACEMENT OF THE EDGE OF THE NEUTRAL SHEET A 35 

HAIL IS THE INCREASE OF THE TAIL FIELD A 36 

IF NO STORM EFFECT IS WANTED, SET ISTQRM=0 A 37 

A 38 

READ (5,40 ISTORM, WIND, STORM, HAIL A 39 

A 40 

IF NU TRACING AT THE OPPOSITE MERIDIAN IS WANTED, SET N0$PLT=1 A 41 

OTHERWISE SET NCSPLT=0 A 42 

A 43 

READ (5,40) NOSPLT A~ 44 

A 45 

GAMMA IS TOTAL ENERGY IN UNITS OF REST MASSES A 46 

A 47 

READ (5,39) GAMMA A 48 

A 49 

FL0NG=1P0. IS NCCN MERIDIAN A 5D 

ALT MUST BE GIVEN IN EARTH RADII FROM CENTER OF FAR tw A 51 

A 52 ' 

READ (5,39) ALT , FL AT , FLONG A 53 

IF (ALT. LT. 0.5) STOP A 54 

KQNT RL=0 A 55 

COS P I T= 1 • A 56 
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KMIN=i A 57 

GSG=GAMMA*GAMMA A 58 

GSCM1=GSQ-1. A 59 

ORLCNG=FLONG+START A 60 

IF I IST0RM.EQ.1 ) 0RL0NG=FL0NG+180. A 61 

A ALT = AL T A 62 

EKMAX = i\MAX a 63 

0CG$= 1 . / EKMAX A 6a 

DO 2 K= 1 , KM AX A 65 

ENERGY! K) =GAMMA-1. A 66 

RLGSTTI K ) = • J . A 67 

2 SDR I FT ( K ) =•’ . A 68 

V ( 1 ) = AL T A 69 

V(i)=<9 J.-FLAT1/57. 2957795 a 7P 

V(3)=FLCNG/57. 2957795 A 71 

0Uf>il = vm A 72 

DUM2= V 1 2 ) A 73 

DUM3= V ( 3 ) A 74 

A 75 

DEFINE INITIAL FIELD LINE A T6 

A 77 

CONTINUE A 78 

CALI. INVAR LV( 1 ) , V ( 2 ) , V( 3) , ERR,BDIM< 1 ) ,FIDIM<1 ) ) A 79 

IF (K.0NTRt_.EC. . > GO TO A A 80 

WRITE (6, All FLONG, FLAT, ALT A 81 

GO TC 5 A 82 

WRITE ( 6 , A 2 ) FLCNG, FLAT, ALT A 83 

DO 6 J= 2 , JUP A 8A 

0LAT = 9:_r.-VN2< J>«57. 2957795 A 85 

GLCNG=VN3< J>»57. 2957795 A 86 

IF (OLGNG.CT.ie:. ) OLONG=OLONG-36C. A 87 

IF (GLONG.LT. ( -130. ) ) 0 L ON G = OL ON G + 3 6 • A ff8 

WRITE (6, A3) J,VNi< J ) , GL AT , OLON G , B < J) A 89 

IF (VNl ( J) .GT.12. ) GO TO i A 90 

CONTINUE A 91 

ERRl=J.‘i*ERR a 92 

IF (KGNTRL.EQ.l) GO TO 29 A 93 

N= A A 9A 

IF (JUP.LE.A) N=3 A 95 

CALL INVAR (VN1(N),VN2(N),VN3(N),ERR,B2,FI2) A 96 

SPATHIl 1 =F 10 IM( l) + (FIDIM(lI-Fl2l*2.*BDIMIl)/(BDIMm-B2) A 97 

call EQUAT (VNEARm,VNEAR(2),VNEAR(3) , EQB ,fc ALT ,ELAT ,E LONG TERM ) A 98 

IF (EALT. GT.12. ) GO TO 1 A 99 

E JALT = C ALT A 100 

VNl 1 - VI'I E AR ( 1 ) A 101 

VN13-VNEARI3) A 102 

RD I PGL = VNl 1 A 103 

0=0.31/ ( RD I P0L««3 . A 1° A 

QG=G/EOB A 105 

SPITCHt 1)=SQRT( l.-EQB/BDIMI ID A 106 

SALT ( 1 ) = ALT A 107 

SLAT l 1 ) =FL AT A 108 

SLONGtl ) =FLONG A 109 

A 110 

LCCP TO DETERMINE ALL I-BM POINTS ON THE INITIAL COMMON FIELD LIN A 111 
CORRESPONDING TO EQUALLY SPACED VALUES OF COSIPITCH ANGLE") AT 112 
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C A 113 

IF UMAX. EG. 1) GC TO 8 A 114 

DO 7 K= 2 , KM AX A 135 

COSPIT=COSPIT-CCCS A 116 

SSQ*l.-COSPIT*CCSPlT A 117 

IF (SPITCHtl ) .LT. COSPIT) SSQ = 1 .-SP I TCH ( 1 )* SP I TCH ( 1 ) A 113 

BO= EQB/ SSQ A 119 

CALL BESECT ( DUM 1 » DUM2 » DUM 3,BALT,BLAT f 6 LONG * ERR 1 ) A 12G 

JJ=JUP-1 A 121 

DUMi-VNK JJ) A 122 

DUM2= VN2 ( J J ) A 123 

DUM3=VN 3 ( J J > A 124 

SALT(K)=VB0(1) A 125 

SLAT (K) *90. -VBO( 2) *57.2957795 A 12o 

SLGNG(K)=VBO(3)*57. 2957795 A 127 

IF (SLONG(K) .GT.180. ) SLONG ( K ) =SLONG ( K ) -360. A 128 

IF (SL0NGIK).LT.(-13C.) ) SLONG ( K > = SLONG ( K ) + 36C. A 129 

CALL INVAR < V BO (1) , VBO < 2 ) ♦ VBO < 3 ) t ERR ,BD I M ( K ) , F I D I M ( K ) ) A 133 

N=4 A 1 3 i 

IF (JUP.LE.4) N = 3 A 132 

CALL INVAR ( VN1 ( N) , VN2( N ) * VN3( N ) «ERR »B2 »F 12) A 133 

SPATH(K)=FIDIM(K)+< F I D I M ( K ) - F 12) *2 . * BD I M ( K ) / ( BD I M ( K ) -B2 ) A 134 

7 SPITCH(K)=SQRT( l.-EQB/BDIM(K ) ) A 135 

8 CONTINUE A 136 

ERR= ERR *40 • A 137 

E R R 1 = ER R *0 • 0 0 L 2 5 A 13b 

ERR2*ERR*C;.025 A 13^ 

A 14 0 

LOOP TO DETERMINE DRIFT VELOCITY ON THE INITIAL LINE A 141 

TO DETERMINE EQUATORIAL DRIFT VELOCITY IN CM/SEC MULTIPLY DRIFT BY A 142 

2.5766 E 05 • (RCIPUL**3) • (MASS IN ELECTRON MASSES) • (GAMMA) • A 143 

• ( BET A*»2 ) . BETA = V/C. A 144 

A 145 

DO 9 K= 1 » KMAX A 146 

F I S T A R = F I D I M ( K ) * C • 9 8 A 1 47 

OALT = SALT ( K ) A 148 

OL AT = SL AT ( K ) A 149 

GLONG^S LONG ( K ) A 150 

CALL SEARCH ( 0 ALT , OL AT rOLONG , BD IM ( K ) , F I STAR , BB , F I , ERR , ERRB , ERR I , RR A 151 

1) A 152 

IF (RK.GT.. 5) GC TU 9 A 153 

CALL EGUAT ( VNE AR ( 1 ) , VNE AR ( 2 ) * VNE AR ( 3 ) , EB ,E ALT , E LAT , E L CNG , E RR1 ) A 154 

DIST = ABS(VNEAR( D-VN11) A 155 

SDRIFT ( K ) = - ( FISTAR-FIDIM(K) ) *QQ / ( D I ST* SPATH ( K ) ) A 156 

9 CONTINUE A 157 

DO 10 K - 1 * KMAX A 158 

IF (SDRIFT(K) .LT.l.E-05) SPATH{K)*0. A 159 

10 CONTINUE A 160 

WRITE (6,44) A 161 

WRITE (6,45) (K#BDIM(K) t FIDIM(K)fSPITCH(K) *SDRIFT(K) tSPATH(K) , EQB , A 162 

1SALT(K) ,SLAT(K) , K = KM IN , KMAX ) A 163 

IF INGSPLT.EG.l ) GO TO 24 A 164 

11 CONTINUE A 165 

FLCNG=URLONG a 166 

12 K A$E= 1 A 167 

DO 13 K= 1 , KMAX A 168 
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RCR I FT ( K ) = ‘ • A 169 

tCALT (K ) = i±l. Hi A i7.) 

GCLCi\GlK) = UIl.il A 171 

RPITCH(K)«0. A 172 

13 eU3H(K)=^. A 173 

A 174 

LOOP TO DETERMINE I-BM POINTS AND SPLIT FIELD LINES AT NEW A 175 

LONGITUDE OR LONG A 176 

A 177 

DO .7 K = K M I N f KMAX A 176 

If (kLuSTT(K) • G T . * i . 5 ) GO TO *6 A 179 

CALL SEARCH I ALT, FLAT, FLG.NGf DDIM(K) ,FIDIM(K) »BB.FI ,E*RfERRB,FRRI ,R A 183 
iLCSTT (K )) A 131 

IF (RLOSTT(K) .GT.v.5) GO TO 16 A 182 

'WRITE (6,46) K , d H I M ( K ) A 1S3 

IF ( { At- S ( FLONG ) • GT . 1 . ) . AND . ( ( i 8 J AB S ( FLUNG ) ) . G T. 1 . . ) ) GO TO 15 A 184 

OG 14 J = 2 j JU P A 155 

UL A T -V' • -V N 2 ( J ) * : 7 .2957795 A i36 

OLCN’G=y N 3 ( J )* 57 . 2^5779? A 137 

IF (LLONG.GT. l c 1 . ) CJL0 nG=GL 0NG- 3b ; • A IFF 

IF (CLCNG.LT . (-IN. . ) ) CLONG = OLONG + 36 1. A 1*9 

WRITt (6,43) JtVNltJ ) , OL A T , 0 L ON G , B ( J ) A 190 

14 CONTINUE A 191 

15 CONTINUE A 192 

N=4 A 193 

IF (JUP.L2.4) N = 3 A 194 

CALL INVAR ( V 4 1 ( N ) , VN'2 ( N ) , VN 3 ( N ) , E RR 2 , B 2 , F I 2 ) A 195 

RPATh(K)=FIOIM(K )+(FiniM(K)-Fi: ) * 2 . *BD I M ( R ) / ( RD l P ( K ) -RE ) A 196 

CALL cWAT ( VNE Ak( 1 ) , VNEARt 2 ) , VNCAR t 3) , C-CB.3 ( K) t EOALT ( K ) , ELAT, ELLON A 197 

l b ( K ) , E R R 1 ) A 198 

HALT (K) *VSAVE ( i ) A 199 

RL AT (K) = FL AT A ? ,N 

RLCNG ( K ) -FLONG A Z";l 

RPITCFt K ) =SQKT( 1 .-CCIMM K ) /Q01M { K ) ) A 2'2 

IF (K.FO.KMAX) GO TO 17 A 233 

flL = RClM (K + l ) A c , <+ 

IF (£QBB(K) .GT.BO) 60= ( BD 1 1‘ ( K ) + FbBD ( K ) ) * ; . 5 A 2 5 

CALL BE SECT ( VS AVE ( I ) , VSAVU £ ) , VSAVt ( 3 ) ,UAL T ,6LAT,PL0NG,£RR1 ) A 2! 6 

ALT * V 50 ( 1 ) A I 'l 

FL AT = BL AT A2CB 

FLCNG=F L GNG A 109 

GO TO 2 7 A 21' 

16 IF (K.lQ.KMIN) KA$E=: a 211 

KK=K A 212 

IF (KASe.EC.l) GC TO 19 A 213 

IF (K.CU.KMAX) CC TU 17 A 214 

K1=K+1 A 215 

ALT=SALT ( K 1 ) - 1 • A 21b 

IF (ALT .LT.l. ) ALT = 1 . A 217 

IF ( ALT .GT . 5. ) A L T = 5 . A i;r 

FLAT = A6S(SLAT(Ki ) )-5. A *.19 

IF ( ABS(SLAT(K1) ).LT.5. ) FLAT=*6. A 22. 

17 CONTINUE A “21 

GO K = KM IN, KM AX A 222 

IF { RLuSTT ( K ) .LT . ‘ . 5 ) GO Tu 21 A 223 

13 CC.\THUF A 2 24 
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o o o 


19 

2G 

21 


C 


23 


25 


GO Tl * A 225 

OG Z K=KK,KMAX A 226 

KLCSTT ( K )= I • A 227 

CONTINUE - A 218 

A 229 

LCOP TG DETERMINE DRIFT VELOCITIES AT SPLIT FIELD LINES A 2 3;? 

A 231 

DO 22 K=KMIN,KMAX A 232 

IF (RPATH(K) .LE.l.E-iO GO TO 22 A 235 

IF (RLUSTT(K) .GT.w.5) CO TO 22 A 234 

VN l = GO ALT (KI A 235 

V N 1 3 = E 0 L C N G ( K ) / 5 7 . 2 9 5 7 7 9 5 A 2 36 

ULi = C/rOHB(K) A 237 

FISTAR=FICIM(K) #2*98 A 230 

UALT = RALT ( K ) A 2.39 

ULAT = i<L AT ( K ) A 24.) 

0 L L N G -H L Li i4 G ( K ) A 241 

CALL SEARCH ( UALT , CL AT , OLONG , BD IM t K ) , F I S TAK , BB , F I , ERR, ERRB » ERRI ,RR A 24^ 

1 ) A 243 

IF (RR.GT. '■ * 5 ) GC TO 22 A 244 

CALL cQUAT ( VNEAR( 1 ) , VNEAR( 2 ), VNEAR t 3) , E6 , E ALT* ELAT ,ELCNG ,ERR1 ) A 245 

G I S T = A 3 S ( V N E A R ( I ) - V N 1 1 ) A 246 

RCRIFT A NO SCRIPT HAVE AS COMMON FACTOR *; • 3 1 / t RDI PGL* * 3 ) A ^47 

KD.R I FT l K ) = — ( FISTAR-FIDIM(K) ) *OG / ( D I $T*KPATH ( K ) I A 2 4R 

CG is T 1 1 vU ?: A 2 49 

CO 2 3 K = KM l N $ KM AX A 2 5 

IF (RCkIPT(K) .LT.l.E- ,5) R P A T H ( K ) = J * A- 251 

CGNT PJOE A 2 52 

rtiUTL lt>,47) AALT. FRONT, SHEET1,SHEET2»FSHEET A 253 

K:<ITL (*5*40 A 254 

WRITE (5, 45) (K, FDIM(K )* FIDIM( K ) ,RPITCH(K) t RDRI FT(K) t RPATH(K) , FOBB A 255 

i (K) *KALT(K) ,RLAT(K) ,K=KMIN,KMAX) A 256 

WRITE ( 6,4Q) A 257 

WRIT?. (5,3 ) A 2 58 

WRITE U,-:n (K,EGLUNG(K) , C OALT,EOALT(K) ,SPITCH(K) ,RPITCH(KJ ,SDRIF A *59 
IT ( K ) , RLR IFT ( K ) , ENERGY I K ) ,K=KM IN,KMAX ) A 2 60 

IF ( I ST Ga m • EG • ) Gu K 28 A 261 

IF (KJNTRL.EO. ) GO To 24 A 262 

Gu TC (25,27,3-), KOmTKL A 263 

FRONT - F R C i si T -W I jC A 264 

SHEET 1=SHE ET1-STGR0 A 265 

FSHEET = FSHEET + H A I L A 266 

KCis T R L - 2 A E 67 

EKr = HKK * . , 2 5 A 263 

wRIT_ h,5:) FRONT T SHEET!, FS MEET A 2^9 

GC TC 7* A 2 7 2 

FRONT -FRONT +V» I\C A 271 

SHEET1 = SHCET 1 + S T C R M a 272 

F$HELT=FSHEET-hA IL a 273 

ORLCNC=CRLCNG+l e O A 2?4 

KG\TKL=E A 2^5 

Cu 2 c K = K M I N , K M A X A 276 

IF { kLCSTT (K).CTO*5) 00 TO 26 A 277 

FISG=FIDIM(K)*riCIM(K) A 278 

CC=ENEH£Y(K)+1. A 279 

CCNST=(CC*CC-1. )«FIS0 A 230 


43 



SOB I = 


BD I M ( K ) * F I SQ 


A 281 

C A 232 

C LOOK St ARCHES FOR THE SHELL HAVING A PREFIXED VALUE OF THE A 283 

C THIRL' INVARIANT, COMPATIBLE WITH CONSERVATION OF THE OTHER TWO A 284 

C A 285 

CALL LOOK (SALT(K),SLAT(K)»SLONG(K)»SQBI,' FLX(K) , BDIK(K) , FIOI A 286 

1M(K) ,ERR,ERRI , FRRB, ERRF, PHI .RLOSTTIKI I A 287 

FI SQ = F I 0 I M ( K ) «F IC IMt K ) A 288 

GAMNEW=SCRT( i . + CCNST/F I SO ) A 289 

ENERGY t K ) =GAMNEW-1 . A 29C 

26 CONTINUE A 291 

WRITE (6,63) FRONT, SHEET 1, FSHEE T A 292 

IF (RLCJSTTI 2) .GT.0.5) GO Tu II A 293 

alt = salt i 2 ) a 2 94 

FLAT = SL AT ( 2 ) A 295 

GO TU 11 A 29fc 

27 CRLCNG=ORLCNG+18C. A 297 

KONT RL= 3 A 298 

IF (NCSPLT.EC.l > GO TO 38 A 299 

GO TL 11 A 300 

23 CONTINUE A 301 

ORLCNG=ORLONG + rSLT A A 3 02 

IF (CKLONG.GT .TERM ) GO TO 38 A 3C3 

FLCNG=ORLONG A 3;4 

IF (RLOSTT(l). GT.0.5) GO TO 1? A 305 

ALT = l<ALT 1 1 ) A 306 

FLAT = RL AT I 1 ) A 307 

GC TU 12 A 308 

C A 309 

C POINTS HAVING THE PREFIXED VALUES OF CONST* ( I •* 2 ) • ( GAMMA**2-1 . ) A 313 

C AND (1**2) *Bf ARE DETERMINED ON THE DISTORTED A 311 

C FIELD LINE A 312 

C A 313 

29 FI — FIOIMt 1) A 3 1.4 

8B=BC[M( 1 ) A 315 

F ISO* F I * F I A 316 

CONST*F IS0*GSQM1 A 317 

COM P = F I SQ*BB A 318 

KM I N = 2 A 319 

K=KM I N A 320 

31 FISC=FIDIM(K)*FIDIM(K) A 321 

CONST=F IS0*GSCM1 A 322 

SQBI=BDIM(K ) *F I SO A 323 

31 IF ( COMP.LTTSQB TT GO TO 36 A 324 

DO 32 J *4, JUP A 325 

IF (BIJ ) .GT . B ( 3 ) ) GO TO 33 A 326 

32 CONTINUE A 327 

33 JUP=J-1 A 328 

IF (JUP.LE.4) STOP 1 A 329 

Si= VN1 ( 2 ) A 3 3D 

S2=VN2( 2 ) A 331 

S3 =VN3 1 2 ) A 332 

S I = F I A 333 

SCOMP=COMP A 33* 

DO 34 J=1,JUP A 33! 

VN1( J)=VN1(J*1) A 336 
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V:M-(J) = VM2tJ + ll A 337 

VnJ ( J)=VN3< J + l) & 338 

Q(J)=B( J + l) A 339 

8LCG(J)=AL0G(B(J) ) A 340 

34 ARC(J)=ABS< ARC! J+l) ) A 341 

JlP = JUP-1. A 342 

DC 33 J = 2 » J C P A 343 

ASUM=AkC(J)+ARC(J+l) A 344 

DX=BLCG(J-1)-8LQGIJ) A 345 

DN=ASUV*AKC(J)»ARC(J+1) A 346 

BCIX (BLOG! J-i)-PLCG( J + l) )«ARC(J)**2-DX*ASUM**2)/DN A 347 

CCL = l DX + AP.CI J + l )-( RLUGI J 1-bLOGt J + l ) )*AP,C ( J ) ) /ON A 348 

SA=.75*ARC( J) A 349 

SC=SA+. 25+ASUM A 350 

0CC=PL0G< J-1)-CCC*SA*SC A 351 

ECCI J)=BCG+CCu»(SA+SC) A 352 

BtG(J)=EXP(0CU+EC0(J)*.5*ARC(J) ) A 353 

35 BEND! J) =EXP< OCC+ECGI J ) • . 5* ( A SUM+AkC < J ) ) ) A 354 

B E C- 1 J U P ) = R E Nl D I J E P ) A 355 

BENCI Jl'P) = B( JUP) A 356 

ECC(JUP) = (2.' / ARC I JUP ) )*ALOG(BEND(JUP)/BEG(JUP) ) A 357 

CALL InTEG ( ABC, BEG, REND, B.JEP, ECO, rLINT) A 358 

FI=FLI -4T A 359 

BB= 8 ( 2 ) A 36" 

FISC=FI*Fl ‘ A 361 

COMP - F I SC* BB A 362 

GO TC 31 A 363 

36 SLCPE=( SCB l-SCOMP )/( COMP-SCGMP ) A 364 

SALT(K)=S1+(VN1(2)-S1)»SL0PE A 365 

&V2(K)=S2+(VN2( 2 ) -S2 ) * SLOPE A 366 

GVi(K)=Si+(Vk3( 2)-S3)*SLOPE A 367 

S I =S I + ( FI-SI ) * SLCP c A 368 

SLAT (K) = 9 «.-C,V 7 (K)*57. 2957795 A 369 

SLCNCIK)=QV3(K)*57. 2957795 A 370 

IF I SLGNG I K ) .GE.18-' . ) SLONGI K ) =SLONG ( K ) -363 . A 371 

IF I SLONGI K).LT.-18-'*.) SLONG I K ) =SLONG l K ) +1 80. A 372 

FISC-SI *S I A 373 

GAMNEW=SGRT( 1 . +CCNST/F I SU ) A 374 

ENERGY I K ) = GAFIN F|«-l , A 375 

FICIMK )=S l A 376 

S I T= ABS I S I N I QV2 I K ) ) ) A 377 

CALL MO DM AG I S ALT l K ) , S l T , GV 3 I K ) , BR , BT , BP ,B6B ,QV2 ( K ) i A 378 

BDIMIK) =BBB A 379 

K=K+'. A 380 

IF (K.LE.KMAX) GC TO 30 A 381 

CALL EtiUAT I VNE AR 1 1 ) , VNE AR ( 2 ) , VNEAR I 3 ) , EQB , E ALT , ELAT ,ELONG , ERR1 ) A 382 

DO 37 K = KM IN, KMAX A 383 

RLOSTT I K ) =v . A 384 

CALL FLUX (SALT, ELAT, ELONG.SALT(K) , SLAT I K) , SLONG I K ) .BDIMIK) .FIDIMI A 385 

IK) ,FLXIK ) , ERR, ERRB, ERR I , RLOSTTI K ) ) A 386 

37 SPITCH(K)=SQRT( l.-EClB/BDIMIK) ) A 387 

WRITE (6,54) A 388 

WRITE (6*45) IK, BDIMIK), FIDIMIK) .SPITCHIK) ,FLX(K) .ENERGYIK) ,EALT,S A 389 

1ALTIK), SLAT(K),K=KMIN,KMAX) A 390 

ERR= ERR *4 . A 391 

CRR1=lRR*0.5w J25 A 392 
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A 393 
~A $94 
A 395 
A 396 
A_ 397 
A '398 
A 399 

tt w 

A 401 
A 402 
A 403 
A 404 
A 405 
A4 TS 
A 407 

44 FORMAT (IX / / / 2X , 1HK , 9X , 4HBD I M , 1 2X , 5HF I D I M , 9X , 6H SP I TCH , 9X i 6HS0R IFT t A 408 


19X,5HSPATH,11X,4FE06 , 9X , 4HS AL T , 4X , 4H SL AT / ) A 409 

45 FORMAT 1 I3.F15.6.4F15.4, F15. 6, 5X, 2F7. 2///) A" 410 

46 FORMAT (/21H SPLIT FIELD LINES/12 0»F2 0.8/1 A_4U 

47 FORMAT ( IX / / 50X , 5 F 10 . 2/ / ) ' M2 

48 FORMAT ( 2X , 1HK, 9X» 4HRDIM, 12X , 5HF IDI M ».9X , 6HRPI TCH »9X ,6HRDR1FT »9X»5H A 413 

IRPATH,11X»4HEQGE»9X,4HRALT»4X, 4HRL AT/ ) A 414 

49 FORMAT ( 28 X , 26HSH E L L SPLITTING PARAMETERS//) A 415 

50 FORMAT ( 2X , 1HK , 4X , 6HEQLONG, 12X » 11HE ARTH RAD I I , 16X , 9HC0S P1'TCH,13X, A 416 

114HDKIFT VELOCITY, 1 3X , 6HENERGY / ) A 417 

51 FORMAT ( I3,Fi:.2,5X,2F10.3,5x,2fl0.3,2El5.3,F14.6//'/T A 418 

•52 FORMAT (2X//77H WIND AND STORM ARE R AG I NG/50X , 3F20 . 2/ ) A 419 

53 FORMAT (2X//21H SPACE IS CALM AGA IN/50X, 3F 20. 2 / 1 A 420 

54 FORMAT (lX///2X,lHK,9X,4HBDIM,12X t 5HFIDIM,9X,6HSPITCH,9X,6H FLUX , A 421 

19X.5HENCRG, 1 1 X , 4HE AL T , 9X , 4HS ALT , 4X , 4HSL A T/ ) A 422 

END A 423- 


ERR2=EKR*0.025 
ALT = S AL T ( 2 ) 

FL AT = SL AT ( 2 ) 

FLCNG*SLONG( 2 ) 

IF (NCSPLT.EQ.l) GO TO 25 
GO TO 11 

38 ERR = ERR*0.:>25 
GO to L 

C 

c 

39 FORMAT (4F13.4) 

40 FORMAT ( I5.3F10.2) 

41 FORMAT ( 29 H DISTORTED COMMON FIELD L I NE /40X , 3F 1 5. 2 / / 1 

42 FORMAT (1X////27H INITIAL COMMON FIELD L I NE /40X , 3F1 5. 2 //I 

43 FORMAT ( 1 5 , F15 . 4 , F 10 . 2, F 1 5 . 2 , E20 . 5 ) 


0 0,00 


SUBROUTINE SEARCH ( ALT ..FLAT, FLONG, SB ,S I . BB , F I , ERR , ERRB , ERRI , RLOST ) B 1 

COMMON B(200 ) ,VNl(l200> ,VN2l2Co) ,VN3(20C ) .ARC120C) ,VNEAR(3) ,VB0(3) , B 2 

1VSAVEI3) ,BO,BNEAR*JUP,'MMM B 3 

COMMON FRONT, SHEET 1 , SHEET2, F SHEET B A 

B 5 

SUBROUTINE DEFINES FIELD LINE GOING THROUGH POINT OF PREFIXED B 6 

B AND SECOND INVARIANT I , AT A GIVEN LONGITUDE ‘ """ " 3 7 

B 8 

DIMENSION V ( 3 ) , V 1 ( 3 ) , V2I3! 9 9 

DV=*.i •02 B ID 

RLOST=j. B II 

MCH£CK=9 B 12 

ICHECK=0 B 13 

SERR=ERR B 14 

SERRB=ERRB B 15 

S ERR I = ERR I B 16 

V ( 1 ) =ALT B 17 

IF (V(l) .LT.C.5) GO TO 11 8 13 

V(2)=(9C.-FLAT)/57. 2957795 B 19 

V(3)=FL0NG/57. 2957795 9 20 

OCLT=i . 5708 B 21 

I CON= 1 B 22 

1 I L I T = 1 B 23 

DELV2=DV B 24 

ICHECK= ICHECK+1 B 25 

IF ( I CHECK .GT • 2\i ) GO TO 12 B 26 

2 S I T = A B S ( S I N ( V (2 ) ) ) B 27 

CALL MUDMAG ( V tl) , S I T , V ( 3 ) , BR, BT , BP , BB , V I 2 ) ) B 23 

3 FAC=1.-(SB-BB)/(3.*SB) B 29 

IF (FAC. GT. 1.5) F AC= 1.5 B 30 

IF (FAC. LT. 0.666) FAC = 4.6e>6 B 31 

V(i)=V( 1)*FAC B 32 

IF <(V<1). GT.KJ.). OR. <V!1). LT.C.5!) GO TO 11 E 33 

VI ( 1 ) =V ( 1 ) B 34 

VI (2 ) =V ( 2 ) B 35 

MCHECK=MCHECK+1 B 36 

IF ( MCHECK . GT . 1 5 ) GO TO 13 B 37 

CALL MODMAG ( V ( 1 ) , S I T , V ( 3 ) , BR, BT , BP , BB , V ( 2 ) ) B 38 

IF ( A BS < (BB-SB)/SB).GT.ERRB) GO TO 3 B 39 

MCHECK= j B 40 

IF (ILIT.NE.l) GO TO 7 B 41 

I L I T=2 B 42 

CALL INVAR ( V 11 ) , V ( 2 ) , V ( 3 ) , ERR , BB , F I ) B 43 

IF (JUP.LT.j) GO TO 14 B 44 

V2 ( 1 ) =V ( 1 ) B 45 

V2 ( 2 ) -V ( 2 ) B 46 

B2=BB B 47 

F 1 2 = F I B 48 

IF l ABS( (FI-SI )/Sl ) .LE.ERRI ) GO TO 8 B 49 

IF (A BS(V(2)- DCLTJ.LT. 0,1) GO TO 6 B 50 

SGN=SIGN( 1., ( FI-SI ) ) B 51 

IF (VI21.LT. DCLT ) GO TO 4 B 52 

DELV2=-SGN*DELV2 8 53 

GO T 0 5 B 54 

4 DELV2=SGN*DELV2 B 55 
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5 

V(2I*VUHDELV2 




B 

56 


GO TO : 




B 

57 

6 

V { 2 ) =V ( 2 ) +CEL V2 




B 

58 


CALL INVAR ( V ( 1 ) , V ( 2 ) , V { 3 ) , ERR , BB , F I ) 




b 

59 


IF (JLIP.LT •'■'?) GC TO 14 




B 

60 


IF {(FI -SI )*( FI-FI2) .LE.O. ) GO TO 2 




B 

61 


V ( 2 ) =V { 2 )-2.*DELV2 




B 

62 


GO TU 2 




B 

63 

7 

HI = BH 




B 

64 


CALL INVAR (V(i),V(2)tV(3).ERR,BB.FI) 




B 

65 


IF (JUP.LT.i-) GC TO i 4 




B 

66 


IF ( AES ( (Fl-S I ) /SI > .LE.ERRI > GO TO 8 




B 

67 


FACT=(S I -F I )/(Fl2-FI ) 




B 

68 


IF ( A0S ( FACT ) .GT.3. ) F AC T= 3 . * S I GN ( 1 . , F AC T ) 




B 

69 


V(l)=Vi ( 1 ) + (V2( 1)“V1( 1) J*FACT 




B 

70 


V (2 >=V1 ( 2) + ( V2( 2)-Vl( 2 > )*FACT 




B 

71 


Y=AMINi ( A8S(V(? )-Vl( 2) )« ABS( V(D-V2(2) ) ) 




B 

72 


IF (Y.GT.ABSt V1(2)-V2U> ) ) GO TO 1 




B 

73 


OV = Y 




B 

74 


GO TC : 




B 

75 

6 

CONTINUE 




B 

76 


IF { ICON. EG. 2) GC TO 9 




B 

77 


I C C N = 2 




B 

70 


DV = L tf*, . 1 




B 

79 


ERK=2RK*v. . 325 




B 

80 


ERRB*ERRB*C.. d 




B 

01 


ERRI=ERRI*' .< 4 




B 

82 


GO TC i 




B 

83 

9 

AL T = V { 1 ) 




R 

84 


FLAT = 9 v • -V ( 2 ) *5 7.295 7795 




B 

85 


FLONG=V( 3)*57. 2957795 




B 

86 


IF (FLUNG. GT. 180.) FLONG = FLONG- 36G . 




B 

87 


IF ( FLONG.LT. (-1*0. > > FlONG-FLGNG+362 • 




B 

80 


CO 1C 1=1,3 




B 

89 

10 

VSAVit I )=V( I ) 




B 

90 


GO TC 16 




B 

91 

11 

WRITE (6,17) V( 1), FLAT, FLONG 




B 

92 


GO TC 15 




B 

93 

12 

WRITE (6,18) ICPECK 




B 

94 


GO TC 15 




B 

95 

13 

WRITE (6,19) I CP ECK , MCHECK 




B 

96 


GO TO 15 




B 

97 

14 

WRITE ( c»2n ) I CHECK, MCHECK, ALT, FLAT, FLUNG 




B 

98 


JUP = 1 




B 

99 

15 

RLCST = 1 . 




B 

100 

16 

cRR=SERR 




B 

1<U 


CRkB=SERRR 




B 

1 02 


ERR I =$F RR I 




B 

1)3 


RETURN 




B 

104 

C 





B 

105 

C 





B 

12-6 

17 

FORMAT (19H ALT OUT OF L IP I TS/ 3F 1C . 3 ) 




B 

107 

18 

FORMAT t 5 1 H SORRY, BUT I CANNOT FIND THAT 

DAMN 

POINT 

IN ICHECK/2 I 1 

B 

103 


1C) 




B 

109 

19 

FORMAT ( 5 1 H SORRY, BUT I CANNOT FIND THAT 

DAMN 

POINT 

IN MCHECK/2 1 1 

B 

110 


13) 




B 

111 

2 ) 

FORMAT (A7H SORRY, BUT POINT IS IN THAT DAMN POCKET/211 CJ»3E15. 5) 

B 

112 


ENU 




B 

113 


48 


n ® >i o vji^u> 


1 


2 


10 

11 


SUBROUTINE LOOK < ALT, FL AT, FLQNG, SQB U TH I R D , B M _ I R* F f M I R , E RR7rR 

1RI .ERRBfERRF, PH I , CUTL ) 

COMMON BC 200 > , VNll 2C0 ) * VN2 ( 200 ) , VN3 1 200 ) f ARC ( 20C ) * VNEAR ( 3 ) * VBO ( 3 ) » 
1 VS AV E ( 3 ) t BQv BNEAR, JUP»MMM 
COMMON FRONT, SHEET It SHEET2, FSHEET 
IF (THIRD. LE. 9. ) GO TO T 
ERRi=ERR*0 .Cl 
C U T L = u • 

I COUNT- 0 
CONTINUE 

CALL SEARCH ( AL T , FL AT , FL ONG, BM I R , F I M I R , BB , F I , ERR , ERR I , ERRB ,RR ) 

IF (RK.GT.0.5) GG TO 7 

CALL EQUAT ( VNE AR , VNE AR ( 2 ) * VNEAR ( 3 ) , EOB , E ALt;, ELAT 7605^0 , WiTlT 

CALL FLUX ( EALT,ELAT,ELONG, ALT ,.FL A T , FL ONG , BB , F I *PHI ,ERR* ERRB , ERRI • 
1CUT ) “ 

IF (CUT. LT. 0.5) GG TU 2 
IF ( I COUNT . EQ .0 ) GO TO 7 
FAC- FAC *0 « 5 
F I M I R=FOLD 
GO TU 5 

DPHI=THIRD-PBI 

IF IABS(DPHI7THIRD).LE.ERRF) GO TO 9 
IF ( I COUNT . GT *0 ) GO TO 3 
FAC=-l./( 12.56* EALT*EQB) 

GO TO 4 

FAC= (FGLD-FIMIR )/lPOLD— PHI ) 

FOLD=FI MIR 
POLC^PH I 
SUM = DPH I *F AC 

IF ( ABS (SUM/FIMIR) .LT.0.1) GO TO 6 

if ( abs(sum/fimir).gt. m . T~rAr= FAC* J7I 
F AC=FAC *C • 5 
GO TO 5 

FIMIR=F IMIR+SUM 
F I S G=F I MIR*FIMIR 

0M I R=SQB I / F I SQ ) 

1CCUNT=IC0UNT+1 

IF ( I COUNT . EG .9 ) GO TO 8 

GO TO 1 

WRITE (6*10 

CUTL= 1 . 

GO TO 9 
WR ITE (6, IT) 

CUTL= 1 • 

RETURN 

FORMAT (24H ‘ NOT ACCESSIBLE IN LOOK) 

FORMAT ( 39H I CANNOT FIND THAT DAMN SHELL IN LOOK) 

END 


C 1 
C 2 
C 3 
C 4 
C" 5 
C 6 
C 7 
C 8 
C 9 


C 10 
C 11 
C 12 


T TS 


C. 14 
C 15 
C 16 
C IT 
C 18 

r it 

C 20 
C 21 
C 22 
X 23 
C 24 

r 2*r 

C 26 
C" 27 
C 28 
C* 29 
C 30 
TT 3T 
C 32 
77 T3 
C 34 
C~ 33 
C 36 

T TT 

C 38 
TT 39 
C 40 
TT 47 
C 42 
T T 3 
C 44 
(T 43 
C 46 
X 47 
C 48 
XT 43 


49 



rn in o 


1 

2 


7 

8 

9 

10 
11 


12 

13 
C 

14 


SU3RCUT INE FLUX { X E ALT ,.X ELAT , XELONG , XAL T , XLAT., XLONG »BM I ft • F I MI R f THI 
IROf ERR, ERRB, ERR 1 , CUT ) 

COMMON B { 20 ) ) * VN1 ( 20 > ) , VN2 ( 200 i »VN3(2C0) , ARC (200) ,VNEAR(3) »VB0(3)« 
1 VS AVE ( 3 ) ,BO.BNEAR f JUPfMMM 
COMMON FRONT, SHE ET1, SHE ET2.F SHEET 
CUT* u 
THIRD*/, . 

KOUNT = u 
CALT=Xt ALT 
E L AT = XL L AT 
E LCNG=X E L ONG 
S ALT*X ALT 
SL AT = XL AT 
S LCNG*XL CNG 
ERR1=ERR*J.31 
0 L * 3 . - • 

CALL VE C p 0T ( E ALT , EL AT , ELONG , APH l ) 

S1=EALT*APHI 
SD=SLCNG/57. 2957795 
IF (SC) 2,3,3 
SO* SC +6 *2331853 07 
GO TC i 

IF (SC-6. 283185307) 5,5,4 

SD*SD-6. 2331853 )1 

GO TC 3 

CONTINUE 

SLCNG=S L GNG + DL 

KOUNT = KOUNT + 1 

IF (KOUNT. GT. 6) GO TO 13 

CALL SEARCH ( S ALT , SL AT , SLONG , 8M I R , F I M I R , BB , F I , ERR , ERRB , ERR I , RR ) 

IF (KR.GT.u.5) GC TU 12 

CALL EQUAT ( ViME AR t 1 ) , VNE AR ( 2 ) , VNEAR ( 3 ) , EQB , E AL T , E L A T , E L ONG , ERR1 ) 
CALL VE C POT ( SALT , EL AT , ELONG , APH I ) 

S2 = E ALT * APH I 
SSD*VNEAR( 3) 

IF ( SSD ) 8,9,9 
SSD=SS0+6. 233185307 
GO TG 7 

IF (SSO-6. 233165327) 11,11,10 

SSO=SSO-6. 233185307 

GO TO 9 

CONTINUE 

DEL* ABS ( SSD-SD) 

SD=S$C 

THIRC*THlRC+( SI +S2 )*0EL*0.5 

si*s: 

GO TC 6 

HKITE (6,14) 

cut*;* 

RETURN 

FORMAT (22H INACCESSIBLE IN FLUX) 

ENC 


0 1 
D 2 
0 3 

D 4 
0 5 

0 6 
D 7 
D 8 
D 9 
0 10 
D 11 
0 12 
D 13 
D 14 
0 15 

D 16 
0 17 

D 18 
0 19 

0 20 
0 21 
D 22 
0 23 

0 24 

D 25 
0 26 
D 27 
D 26 
0 29 

D 3) 

D 31 
0 32 

0 33 

0 34 

D 35 
D 36 
0 37 

D 30 
0 39 

0 40 

0 41 

0 42 

0 43 

D 44 
0 45 

0 46 

D 47 
0 48 

D 49 
0 50 

0 51 

D 52 
0 53- 
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SU8RCUT I NE VECPCT l RO , EL AT , ELONG , A ) N 1 

COHMCN B t 2v0 ) * VNl 1 2CC ) » VN2 ( 20C ) , VN3 I 2dC ) . ARC t 200 ) , VNE AR ( 3 ) ,VB0(3) • N 2 

lVSAVE(3)«BO v BNEAR>JUP«MMM N 3 

CONMON FRONT t SHEET It S HE £ T2 1 F SHEET N 4 

A=u. 31 N 5 

0NG=EL0NG/57. 2957795 N 6 

1 IF (CNG) 2,3,3 N 7 

2 0NG=CNG+6. 283185307 N 8 

GO TC 1 N 9 

IF ( ONG -6*293185307 ) 5,5,4 N 10 

0NG=CNG-6. 283185307 N 11 

GO TO 3 N 12 

5 CGNTINUE N 13 

R= 1 • N 14 

UV={R0-1.)**).G7 5 N 15 

R2 = R + 0V *C * 5 N 16 

VLAT=(9G.-ELAT)/57. 2957795 N 17 

S1T=ABS(SIN(VLAT)) N 18 

6 CALL MOCK AG i R2 , S I T , ONG , BR , BT , BP , BB , VL AT ) N 19 

A= A-ABS ( BT )*R2*CV N 20 

R2=R2+DV N 21 

IF (R2.GT.R0) GC TO 7 N 22 

GO TC 6 N 23 

7 A=A/RO N 24 

RETURN N 25 

END N 26- 
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SUBROUTINE EQUAT ( CUM1, DUMZ, DUM3»EB ,EALT, ELAT, ELONG .ERR) 

COMMON B ( 200 ) , VN 1 C 20 2 ) *VN2 ( 200 J , VN 3 ( 200 ) » ARC ( 200 ),VNEAR(3), VBO ( 3 I . 
1VSAVEI3 ) t BOiBNEAR, JUP,MMM 
COMMON FRONT, SHEET 1 , SHEET2, F SHEET 

SUBROUTINE TRACES FIELO LINE FROM A GIVEN POINT TO MINIMUM B 
MINIMUM B POCKETS AT HIGH LATITUDES ON NOON SIDE ARE IGNORED 

DIMENSION V ( 3 . 3 ) . VN(3), VP(3), Rl<3), R2(3), R3<3) 

ERR1=ERR 
MMM = 1 

jup=; 

V ( 1 * 2 ) = DUM 1 
V ( 2 i 2 ) = DUM2 
V ( i , 2 ) = DUM 3 
ARC C 1 ) = C . 

DCLT=1. 570ft 

ARCU ) = V( 1,2 )*SCRT(ERR)*0.3 
IF ( V ( 2 9 2 ) -DCLT ) 2,3,3 
ARC ( 2 ) =-ARC ( 2 ) 

CALL START ( R 1 , R2 , R 3 . B , ARC , ERR , V ) 

IF tJUP.LT.T) GO TO 5 
DC 4 1=1,3 . 

vp( n=vt 1 , 2 ) 
vn( n=v n , 3) 

CALL LINES < R 1 , R2 t K3 , B . ARC , ERR , J , VP . VN ) 

IF (J.LT.20J) GO TO 7 
ERR=4 • * ERR 
GO T Lj t 
JUP = J 

ERR=tRR*t . 1 

CONTINUE 

WRITE (6, ft) ERR 

GO TU 1 

ERR=hRR 1 

EB=BNEAR 

EALT=VNEARll) 

E L A T = 9u • -V N E A R ( 2 )* 57 . 29 57795 
ELONG=VNEAR( 3)*57.2957795 
IF (EL0NG.GT.16C.) ELONG -ELONG- 360 • 

IF (ELUNG.LT. I-IBC.J) ELONG= ELONG + 360 • 

RETURN 


FORMAT 1 24H ERROR CHANGED IN EQUAT.E15.41 
END 


6 1 
E 2 
E 3 
E 4 
E 5 
E 6 
E 7 
E 8 
E 9 
E 10 
E 11 
E 12 
E 13 
E 14 
E 15 
E 16 
E 17 


E 20 
E 11 
E 22 
E 23 
E 24 
E 25 
E 26 
E 27 
E 20 
E 29 
E 30 

e n 

E 32 
E 31 
E 34 
E 31 
E 36 
E 37 
E 38 
E 19 
E 40 
E 41 
E 42 
E 43 
E 44 
E 41 
E 46- 
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SUBROUT I NE BESECT ( RUM 1 , KUM2 t RUM3 1 BALT , BLAT , BLONG t ERR) F 1 

CGMMUN 3(2.: }) * VNim ) ,VN2(23C ) ,VN3( 200) ,ARC(2‘ C) ,VNEAR{3) , V60(3) , F 2 

lVSAVE(2i ) tBOf BNEARt JUPt MMM F 3 

CCMMCN FRONT, SHEET1, SHEETS, F SHEET F 4 

F 5 

SUBROUTINE TRACES FIELD LINE UPWARDS FROM A GlVEu POINT UNTIL F 6 

A PREFIXED B-VALUE IS REACHED F 7 

F 8 

DIMENSION V<3,2), V N ( 3 ) , VP(3), Rl(3), R2<3), R3<3> F 9 

ERR 1 = ERR F 10 

M M M = 3 F 11 

JUP=_ F 12 

V ( I * 2 ) = RUM! F 13 

V(2,i)=RUM2 F 14 

V(3 f< :) = RUM3 F 15 

A R C ( 2 ) = / • F 16 

DCL T = 1 • 57 j 8 F 17 

ARC(2)=V( 1,2)*SCRT( ERR )*0.3 F 18 

IF { V ( 2 , 2 ) -DCLT > 3,3,2 F 19 

ARC (*.) = - ARC ( 2 ) F 20 

CALL STAKT I R1,R2»R3, B, ARC, ERR, V ) F 21 

IF { JUP . LT . 3 ) GO TO 5 F 22 

DO 4 1=1,3 . F 23 

VP { I ) =V ( 1,2) F 24 

VM I)=V( 1,3) F 25 

CALL LINES ( R 1 , R 2 , R 3 , B ,-ARC , ERR ,.J , VP , VN ) F 26 

IF (J.LT.23?) GO TO "/ F 27 

ERR=4.*ERR F ?R 

GO TO 6 F 29 

JUP=: F 30 

EKR=ERk* 0. i F 31 

CONTINUE F 37 

WRITE (6,8) ERK F 33 

GO TO 1 F 34 

ERR=ERkl F 35 

JUP=J F 36 

BALT=VB0( 1 )*6371.2 F 37 

BLAT=90.-VB0( 2)«57. 2957795 F 3b 

BLCNG=V BG ( 3)*57. 2957795 F 39 

IF ( BLONG. GT. 130. ) BL0NG = BLONG-360 • F 40 

IF (6LOiMG.LT. (-180. ) ) BL0NG= BL0NG + 36 2 . F 41 

RETURN F 42 

F 43 
F 44 

FORMAT (24H ERROR CHANGED IN BE SEC T , E 1 5 . 4) F 45 

ENC F 46- 
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. SUBROUTINE INVAR 1 CUM 1 , DUM2 , DUM 3 , ERR * BB , F I ) 

COMMON P ( 2 ■;«*>) ,VNU2*‘ !■ ), VN2! 2C0 ) , VN 3 ( 200 ) f ARC ( 2 0C ) ,VNEAR<3) ,VB0<3 ) , 
1 VS A Vc ( 3 ) t 80, BNEAR, JUP,MMM 
COMMON FRONT, SHEET 1, SHE ET2, F SHEET 

SUBROUTINES I NV AR , ST ART , L INE S AND INTEG ARE BASED CN MCILWAINS 
INVAR CODE 

DIMENSION V ( ? , 2 ) , V N ( 3 ) , VP(3), BEG ( 20 j ) * BENDI20U), BLOG12JC), EC 
10CN‘), K i ( 3 ) * RH3), l<3(3) 

MMM = .. 

JU° = 1 

Vt 1 ,Z)=CUM1 
V ! 2 » 2 ) = DUM2 
VI j,.) = nUM3 
CRRi=EkR 
ARC ( 1)= . 

ARC(2 ) = V( l,2)*SCRT(ERR)*').3 
DCLT=1.57v3 

IF (V(2,2)-DCLT ) 2,3,3 
ARC ( 2 ) = - ARC ( 2 ) 

CALL START ( R 1 , R 2 , R 3 , B , ARC , E RR , V ) 

IF ( JUP.LT. )) GO TO 8 
DO 4 1=1,3 
VP! I ) =V ( 1,2) 

VN( I)=v ( I , 3) 

CALL LINES ( R i , R2 , R 3 , B , ARC , ERR ,. J , VP , VN ) 

IF ( J .L T • 2 v ' ) GC TO 5 
ERR= ERR *4 • 

WRITE (6,9) ERR 
GO TO ’ 

EH R = ERR 1 
JUP = J 

DO t J=1,JUP 
ARC! J)=ABS(ARC( J) ) 

BLCG ( J ) = ALOG { B ( J ) ) 

JEP= JUP-1 

DO 7 J= Z » J EP 

ASUM= ARC ( J ) + ARC ( J + 1 ) 

DX=BL0G< J-D-BLCGU) 

DN = ASUM* ARC ! J )*ARC( J+l ) 

BCO= ( (BLOG! J-l) -BLOG! J+l) )*ARC( J )«*2-DX*ASUM**2 ) /DN 
CCC=(CX*ARC( J + l ) - C BLOG ( J )-BLOG( J + l) ) * ARC { J ) ) /ON 
S A= * 75* ARC ( J ) 

SC=SA+. 25*ASUM 
DCC=BLUG( J-l )-CCC*SA*SC 
ECU! J)=BCO+CCU*( SA+SC) 

BEG! J) = EXP(DCO + ECO(J )*. 5* ARC ( J ) ) 

BE NO! J) =FXP! DCO+ECC! J )*.5*( ASUM+ARC! J) ) ) 

B F G ( J U P )=BENO(JEP) 

BEND! JUP ) =B( JUP ) 

ECO! JUP)= ( 2.1/ARCIJUP ) )* ALOG ( BEND UUP ) /BEG! JUP) ) 

CALL INTcG ! ARC, BEG, BEND,B,JEP, ECO, FLINT) 

F I = F L I N T 
BB= B ( Z ) 

cont inue 

RETURN 


FORMAT ( 26 H ERROR INCREASED IN INVAR,E15.4) 
END 


G 1 
G 2 
G 3 
G 

G 5 
G 6 
G 7 
G fi 
G 9 
G IS 
G 11 
G 12 
G 13 
G 14 
G 15 
G 16 
G 17 
G 18 
G 19 
G 2D 
G 21 
G 22 
G 23 
G 24 
G 25 
G 26 
G 27 
G 28 
G 29 
G 30 
G 31 
G 32 
G 3 3 
G 34 
G 35 
G 36 
C 37 
G 33 
G 39 
G 40 
G 41 
G 42 
G 43 
G 44 
G 45 
G 46 
G 47 
G 48 
G 49 
G 50 
G 51 
G 52 
G 5 3 
G 54 
G 55 
G 56 
G 57 
G 58 
G 59 
G 50 
G 61- 
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SUBROUTINE START ( K 1 f R2 , R 3 * B , ARC . ERR » V ) H 1 

COMMON B(2 -0 ) , VN1120 * ) , VN2 ( 2 CG ) , VN3 ( 20U ) , ARC (200) * VNEAR ( 3 ) , VB0(3) , H 2 

1 V S A V E ( 5 ) , B0» BNEAR, JUP*MKM H 3 

COMMON FRONT, SHEET!, SHEET2, FSHEET H 4 

DIMENSION 7(3*3), Kl(3), R2(3), K 3 ( 3 ) H 5 

L0CP=1 H 6 

SIT=ABS(SIN(V(2*2 ) ) ) H 7 

1 I F ( V ( 3 , : ) ) 2 , 3 , ? H ? 

2 7(0, 2)=V(3,2)+6. 283185307 H 9 

00 TO : H 10 

3 CALL MO DMA G ( V ( 1 , 2 ) , S I T , V ( 3, I ) ,.d* * 6T , BP , B ( 2 ) , V ( 2 , 2 }) H 11 

R2( l)=BR/6( 7) H 12 

DN= E ( 2 ) *V ( 1, 2 ) H 13 

R 2 ( 2 ) = h T / D N H 14 

R2(3 ) =B P / ( ON *S l T ) H 15 

I S = H 16 

4 DU 5 1=1,3 H 17 

5 V( I,i)-V( I»2)-ARC12)*R2( I) H 18 

SIT*ABS(SIN(V(2 f IP . H 19 

6 CALL MCJDMAG ( 7 , S I T , V ( 3 , 1 ) , BR , BT , BP , B ( 1 ) , V ( 2 , 1 ) ) H 20 

IF ( B ( 1 J - B ( 2 ) ) 7,8,8 H 21 

7 ARC ( 2 ) =-ARC ( i ) H 22 

IF (LCUP.EG.2) GO TO 12 H 23 

LOC P = 2 H 24 

0 0 T (J 4 r. 2 5 

8 R1(1)=BR/B(1) H 76 

ARC ( 3 ) = ARC ( 2 ) h 27 

DN=B(1)*V( 1, U H 28 

Rl(2)=bT/DN H 29 

Ri(3)=BP/(DN*SIT) H 30 

DO 9 1=1,3 R 31 

9 V(I,i) = V(I,2)-ARC(2)*(Rl(I)+R2U)>/2. H 32 

SIT=ABS(SIN(V(2,1) )) H 33 

IS=IS+1 H 34 

SO TO (6,10), IS H 35 

13 30*11=1,3 H 36 

11 f ( I , 3 ) = V ( I , 2 ) + ARC ( 3 ) * { ( 1 • 5 ) *R2 ( I )-.5*Rl( I ) ) H 37 

SO Tt 13 H 38 

12 JUP=-i H 39 

13 CONTINUE H 40 

RETURN H 41 

END “ 42- 
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SUBROUTINE LINES t R 1 , R2 , R 3 , B , ARC , ERR , J , VP , VN ) I 1 

COMMON P(2j‘l),VNHZO . > , VN2 < 2C;.' > , VN3 ( 20. > , ARC 1 2C C > , VNEAR (3) , VBDI 3 > , I 2 
1VSAVEI 3 ) , BO, BNEAR, JUP.MMM l 3 

COMMON FRONT »SHEET1, SHEET2»F SHEET I A 

INTEGER FL AG1 , FLAG 2 I 5 

DIMENSION K 1 ( 3 ) t R 2 ( j ) » R3(3), VN(3>, VP ( 3 ) , R A < 3 ) I 6 

M=MMM l 7 

FLAG!®.' I 8 

DEL = J.i i I 9 

CR0=.._5 l 10 

IF (ERR-, ,.15628) 1,2,* I 11 

I CRE=(ERR**-‘ . 333333333) I 12 

Z A3= ARC( 2 ) I 13 

AAB=AHS(A3) I \A 

SNA= A 3/ A AB I 15 

A. = ARC ( 1 ) I 16 

A2 - ARC ( 2 ) I 17 

AC6 = A3* A3/6... I 18 

VN1 ( 2 ) = V P ( 1 ) 1 19 

VN2 ( i ) = VP( 2 ) I 2,.' 

VN3(2)=VR(3) I 21 

J = 3 I 22 

I L P = 1 I 23 

I S = 1 I 2A 

GO TO 2 I 25 

3 IS=1 I 26 

J = J + 1 I 27 

AC6=A3*A3/5. I 23 

ARC J = Al + A2 + A3 I 29 

AD-IASUM + Al ) / AA I 39 

BO=A$UM./BB l 31 

CD= A i/C C I 32 

A 00 7 1=1,3 I 33 

00=Ri ( I 1/AA-R2I I ) /8B + R 3 ( I )/CC I 3A 

GO TC (5,6) , IS I 35 

5 KT=Ri( I )-( AD» R 1 ( I )-BD*R2( I ) +CD*R 3 ( I )-DD*ARCJ )*ARCJ I 36 

RA(I)=R1(I) I 37 

R1 ( I ) =K2 ( I ) I 38 

R2 ( I ) =R 3 ( I ) I 39 

R3(I)=RT I A3 

V P ( I ) = V N ( I ) I A1 

RBAR=(R2 ( I )+R3I I ) 1/2 .-DO* AO 6 I A2 

VN( I )=VP( I ) + A3«RRAR I A3 

8 I F (VN( 2 ) ) 9, 1 ‘ , U I AA 

9 VM 2 ) =-VN ( 2 ) I A5 

1C IF (\)N(2)-3.1A*692653) 12,12,11 I A6 

II V M 2 )=fc.2831S539?-VN( 2 ) I A7 

GO TL 1 3 I AR 

12 IF ( V N ( 3 1 ) 12,1A,1A I A9 

13 VMj|=VN (31+6.2831353 7 I 59 

GO TO 12 I 51 

1A IF( VN(3)-6. 28 U?53„' 7) 16,16,15 I 52 

15 V N ( 3)=VN( 31-6.2831653 '.'7 I 53 

oO TC 1 A I 5 A 

10 GO TO ( 17, IP' 'S I 55 
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* 


17 


13 

19 


20 

21 


22 

23 


24 


25 

26 
27 


S I T = A BS ( S I N ( VN ( 2) ) ) 

PK E 1 = VN ( 1 > 

PRE*: = PRE1*VN t 2 ) 

PkE3=PREl*SIT *VN ( 3 ) 

CALL MODMAG ( VN , S I T , VN { 3 ) , BR , BT , BP , B { J ) , VN ( 2 ) ) 

R3 ( 1 )=BR/B( J ) 

DN=B( J) *VN( 1 ) 

kj { 2 ) =BT / DM 

K3{^)=RP/trN*sm 

ASUF=A3+A2 

A A = A S U M * A 2 

BB=A3*A2 

CC=ASUh*A3 

I S = 2 

GO TO 4 

3 IT = AB3 (S INI VN( 2 )>) 

IF ( VN ( 1 ) .GT.8. ) GO TO 19 
B ( J ) = B ( J)*( (PRFi/VNt 1) )**3) 

IF (K.EC.I) GO TC 23 

QRT=.5*ACS(R3( 1) )/( . i +ABS ( R 3 ( 2 ) * VN < 1 ) ) ) 

X = ( A B S t V N ( n-PREl)+GRT*ABS(VN(l)*VN(2)-PRE2)+ABS{VMl)*SIT*VN(3)-P 
1RE3 ) ) /( AAB*ERR*SCRT ( l.+QRT*GRT) ) 

GO TC (23,21,23), ILP 

IF (X-3.3) 23,21,21 

A3 = A5*i • 2 * ( 8 • 1 + X )/('.* 8 + X ) 

J = J-i 
I LP= 3 

ASUN=A2+A1 
AA=ASUK*A1 
B B = A 2 * A 1 
CC=ASUr**A2 
CO 22 1=1,3 
V N ( I )=VP( I ) 

R 3 ( I ) = R 2 ( I ) 

R2 ( I ) =R 1 ( I ) 

R 1 ( I ) = k A ( I ) 

GO TU 35 
V N 1 ( J ) = V N ( 1 ) 

VN2(J)-VN(2) 

VN3 ( J ) = V N ( 3 ) 

IF (M.EC.3J GO TC 27 

IF ( ti( J-l) .GT.Pt J) ) GO TO 29 

IF (ABS(VN2( JJ-1.57) .GT.0.26) GO TO 29 

IF (M.EQ.i) GO TC 25 

IF (FLAG1.EQ.1) GO TO 29 

FL AG 1 = 1 

CO 24 1=1,3 

VNEAKl I ) = V N ( I ) 

BNEAR = BI J ) 

GO TC 29 
BNEAR=Bt J-l ) 

DO Z c I=i,3 
VNEARt I ) = V P ( I ) 

GO TO 37 

IF (BCJ)-GT.PO) GO TO 29 
FAC=(BG-B{ J-l) )/(B( J >-R{ J-l) ) 


I 56 
I 57 
I 58 
I 59 
I 62 
I 61 
I 62 
I 63 
I 64 
I ‘ 65 
I 66 
I 67 
I 68 
I 69 
I 70 
I 71 
I 72 
I 73 
I 74 
I 75 
I 76 
I 77 
I 78 
I 79 
I 80 
I 61 
I 82 
I 83 
I 84 
I 35 
I 86 
I 87 
I 88 
I 89 
I 90 
I 91 
I 92 
I 93 
I 94 
1 95 
I 96 
I 97 
I 98 
I 99 

I too 

I 101 
I 102 
I 103 
I 1C4 
I 105 
I 106 
I 107 
I 108 
I 109 
I 110 
I 111 
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CO 26 1 = :, 3 I 112 

28 VBCI I > = V P ( I )+<VN( I )-VP(I) )*FAC I 113 

ARC(J)=ARC( J ) *F AC I U* 

B { J ) = DO I 115 

VNl(J)-VEOll) I 116 

VNt ( J ) - V 80 ( 2 ) I 117 

vn? ( j ) =vbo 1 3 ) i ue 

GO TO 37 I 119 

29 IF IJ.GC.20j) f,C TO 37 I 12'' 

Ai = A2 I 121 

IF (M.NE.v) GO TO 2: I 12 2 

IF { R ( J ) -6 ( 2 ) ) 31 » 2*- 1 36 I 123 

30 I L P-2 I 1 2 A 

A2=A3 I 125 

IF (M.EG.l ) GO TO 35 I 126 

A3=A3*.2»(6.+X)/( .?+X) I 127 

AM=U.-R3(2)*VN(1> )*VN( 1)*CKE I 128 

IF ( ABS ( A3 ) - AM ) 32,32, 31 I 129 

31 A3=SNA»AM I 130 

32 IF <SNA*R3 ( 1 ) +. 5 ) 33,33,35 I 131 

33 AM=-«5»SNA*VNl 1 I/R3I 1 ) I 132 

IF (ABS l A3 )-AM) 35,35,34 I 133 

34 A3=SNA* AM I 134 

35 ARCIJ+l ) = A 3 I 135 

AAB=ABS(A3> I 136 

GO T l. 3 1137 

36 CONTINUE I 138 

37 RETURN l 139 

END I 1*0“ 
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4 

5 

6 


7 

8 

9 

10 

11 

12 

13 

14 

15 


SUBROUTINE INTEG 1 ARE , BEG * BEND, B , JEP ,ECO ,F I ) 

COMMON B(209> *VNIC200) , VN2( 2001 , VN3(200) , ARC (200 ,VNEAR(3) ,VB0(3) . 
1 VS AVE ( 3 ) ,B0, BNEAR, JUP,MMM 
COMMON FRONT, SHEET 1, SHEET2, FSHEET 
DIMENSION BEG ( 200 ) * BEND(200)» EC0I200) 

KK= JEP 

IF IKK— 4) 3,2,4 
KK-KK-1 

A=B(KK-1 I/BI2 ) 

X2 = B ( KK ) / B ( 2 ) 

X3 = B ( KK + 1 ) /B ( 2) 

ASUM=ARC(KK)+ARC(KK+1) 

DN=ARCl KKI *ARC(KK+1)*ASUM 

BB= {-A*ARC(KK + 1 ) *IARC( KK J+ASUM ) +X2* A SUM**2- X3* ARC (KK)**2) /ON 
C= ( A * ARC ( KK 4- 1 )-X2*ASUM+X3*ARC(KK) ) /D N 
F 1=1. 57079 6326* 11 .-A+BB*BB/< 4.*C ) ) / SQRT ( AB S ( C ) ) 

RETURN 

T = SORT ( 1 • -BEND ( 2 ) / B ( 2 ) ) 

FI=(2.*T-AL0G( ( l.+T)/(l.-T)) )/ECU(2) 

IF IBU)-BEND(KK) ) 6,6,5 
KK = KK + I 

T= S CRT { ABSI 1.0-BEG(KK)/B(2) ) ) 

FI=F 1 -( 2.*T-ALGG( t 1 . + T ) / ( 1 T I ) ) / ECO ( KK > 

KK=KK-1 

DO 15 1=3, KK 

ARG1 = 1 * “BEND ( I ) / P C 2 ) 

IF ( ARG1 ) 7,7,8 
TE= 1 • E-5 
GO TO .9 
TE=SQRT ( ARG1 ) 

ARG1= 1 • -BEG( I ) / B { 2 ) 

IF ( ARG1 ) 11,11,10 
T B = SQRT ( ARG 1 ) 

GO TC 12 
T B= 1 • E-5 

IF ( ABStECOt I ) I-2.E-5) 13*13.14 
F I - F I ♦( (TE + TB)*(ARC( I ) +ARC ( I +1 ) ) 1/4. 

GO TU 15 

F I = F I + ( 2.*(TE-T8 ) - ALQGt < 1 . +T E ) * 1 1 TB ) / ( ( I . -TE > * ( 1 . + TB > ) ) ) /ECOU ) 

CONTINUE 

RETURN 

END 
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J 21 
J 22 
J 23 
J 24 
J 23 
J 26 
J 27 
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J 30 
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J 34 
J 35 
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J 40 
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SUBROUTINE MODMAG <RR,5INTH,PPHI,BR.BThETA.BPHI.BB.THET) 

SUBROUTINE MODMAG ( RRi* S I NTH. PPH I , BR, BTHETA , BPH I * BB. THET ) J 1 

COMMON FRONT, SHEET1. SHEET2, 'SHEET J 2 

J 3 

SUBROUTINE ASSEMBLES MAGNETIC FIELD FROM TAIL. MAGNETOPAUSE AND J « 

INTERNAL DIPOLE - TOT ENQUIRIES WRITE To Gil BERT MEAD. GODDARD J 5 

SPACE FLIGHT CENTER, GREENBElT MARYLAND 20771 J 6 

J 7 

DIMENSION GG<7,7> J 8 

DIMENSION G ( 7 * 7 ) , C0NST(7.7>, P<7.7)» DP < 7 , 7 ) , Sp<7), CP<7) J 9 

C0STH«C0S(THET) J 10 

S!NPHI»-SIN(PPWn J 11 

C0SPHI«-C0S<PPHI ) J 12 

ROsFRONT J 13 

RUSHEET1 J 14 

R2sSHEET2 J 15 

BCS«FSHEET J 16 

IF (JFIRST-13) 1,5,1 J 17 

JFIRST.13 J 18 

J 19 

SET UP INITIAL CONSTANTS THE FIRST TIME AROUND J 20 

J 21 

DO 2 N» 1 » 7 J 22 

DO 2 M«l , 7 J 23 

GG(N,M)sO. J 24 

NMAX*7 J 25 

J 26 

THE following COEFFICIENTS ARE STHMIDT-NORMALIZED J 27 

J 2ft 

GG C 2, 1 ) ■ -0 , 2511 1E5 J 29 

GG(3,2)»-0.12424E5 J 30 

GG<4,l)i-0 , 00716E5 J 31 

GG(4,3)>-0,02333E5 J 32 

GG(5, 2)*-0, 02397E5 J 33 

GG(5, 4>i-0, 00163E5 J 34 

GTJC6 » 1 ) «0 , 0O569E5 J 35 

GG(6, 3)i-0. 01078E5 J 36 

GG(6, 5)i-0, 00103E5 J 37 

GQ(7,2)iO .00126E5 J 38 

GG(7, 4)i-0. 00187E5 J 39 

GG(7, 6>i-0, 00041E5 J 40 

P * 1. 1 ) * 1 . 0 J 41 

DP ( 1 , 1 ) i0 . J 42 

SP(1)»0. J 43 

CP(1>«1. J 44 

DO 3 N«3.NMAX J 45 

FN=N J 46 

N2*N-2 J 47 

DO 3 Mil , N2 J 48 

FMeM J 49 

3 C0NST(N.M)i( (F n- 2. )**2-(FM-l, )**2)/( <2.*FN-3. )*<2.*TN-5. ) ) J 50 

DIMENSION SHMIDT (7,7) J 51 

SHMIDT(1,1)«i.o J 52 

00 4 N«2,7 J 53 

FNiN J 54 

SHMIDT <Ni 1 )»SMMIDT<N-1, 1)*< * N*FN-3 . 0 >/ CFN-1 . 0 ) J 55 

FACT i2 , 0 J 56 

DO 4 M»2,N J 57 

FM«M J 58 

SHMIDT(N,M)«SMMIDT(N.M-1)*S3RT( <FN-FM*1.0)*FACT/(FN*FM-2.0) ) J 59 
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4 FACTsl.O J 60 

5 IF (R3-ROOLD) 6.9.6 J 61 

6 ROOLDsRO J 62 

DIMENSION F AC ( 7 ) J 63 

FAC(2)sR0**3 J 64 

DO 7 Nx3,NMAX J 65 

7 FAC(N)sRO*FAC<N-1) J 66 

DO 8 Nx?,NMAX J 67 

DO 8 M*i,N J 68 

9 G(N,M)sSHMlDT(N.M)*G3(N.K)/ r AC(N) J 69 

9 CONTINUE J 70 

0 BEGIN CALCULATION FOR SPEC I r I ED INPUT J 71 

: j 72 

CTsCOSTw J 73 

STsSINTH J 74 

SP(2)*S!NPHI J 75 

C p (2)sC0SPHl J 76 

0 J 77 

0 CALCULATE S I N * M*PH I ) AND COS(M»PHI) J 78 

: j 79 

DO 10 M*3.NMAX j 80 

sp(m)»sp<?)*cp(m-d*:p(2)*s :, (m-i) j ei 

10 C p (M)sCP(2)*CP(M-l)-5P(2)*S 3 (M-l) J 82 

Rsl’R J 83 

A*l, J 84 

RAsR/A J 85 

ROAxl, J 86 

BR=0. J 87 

BTHETAxO. J 88 

B p HlxO. J 89 

FNst, J 90 

C J 91 

C CALCULATE SPHERICAL HARMONICS FOR CAVITY FIELD J 9? 

C J 93 

DO 16 N*2.NMAX J 94 

SUMRsO. J 95 

SUMTxO, J 96 

SUMPsO. J 97 

FM=0, J 98 

C J 99 

C DEVELOP LEGENDRE FUNCTIONS AnD THEIR DERIVATIVES BY RECURSION FORM J 100 

C J 101 

DO 15 Mxl.N J 102 

IF (N-M-i) 13.12,11 J 103 

11 P(N#M)scT*P(N-l.M)-CONST(N,H)«P(N-2,M) J 104 

D p (N,M)*CT*DP<N-l,M)-ST*P(N-l.M)-C0NST(N,M)*DP<N-2,M) J 105 

GO TO 14 J 106 

12 P(N»M)xCT*P(N-l.M) J 107 

D p (N,M)sCT*DP<N-1.M)-ST*P(N-1.M) J 108 

GO TO 14 J 109 

13 P(N,N)sST*P(N-1.N-1) J 110 

D p (N,N)sST*DP(N-1.N-1)*CT*P(n-1.N-1) J 111 

14 CONTINUE J 112 

TSsG(N,M)*CP(M) J 113 

SUMR=SUMR*P(N.M)*TS J 114 

SUMT=SUMT*DP(N,M)*TS J 115 

SUMP=SUmP*FM*P(N,M)*3(N.M)*5p(M) J 116 
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FM=FM+1 . 

J 

117 

15 

CONTINUE 

J 

118 


BR=0R-ROA*FN*SUMR 

J 

119 


BTHETAsBTMETA-ROA*SJ'1T 

J 

120 


B p HI=9PHI*R0a*SUmP 

J 

121 


RQA*R0A*RA 

J 

122 


FNiFN*1 . 

J 

123 

15 

CONTINUE 

J 

124 


bphi=bpri/st 

J 

125 

C 


J 

126 

C 

CALCULATE TAIL FIELD 

J 

127 

c 


J 

128 


RCT=R»CT 

J 

129 


RCT2sRCT**2 

J 

130 


RSC*R*ST*CP<2> 

J 

131 


TOP=R2*RSC 

J 

132 


BOT«Rl*RSC 

J 

133 


BX*-BCS*( ATAN(TOP/RCT) -AT AN (BOT/RCT) )/3. 14159265 

J 

134 


BPHI*BPMI*BX*SP(2) 

J 

135 


BRHO*«BX*CP<2> 

J 

136 


BT=BCS*aLOG( (RCT2*T0 3 *«2)/( RCT2*B0T**2 )) /6 . 28310531 

J 

137 


BRsBR*BRHO*ST-BT*CT 

J 

138 


BTHETAsBTWETA*BRHO*CT*BY*ST 

J 

139 



J 

140 

V 

ADD DIPOLE FIELD To CAVITY r I ElD 

J 

141 



J 

142 


R3=fi**3 

J 

143 


BR=BR-6?000.*COSTH/R3 

J 

144 


BTHETAsRTHETa-31000 . *SINTH/R3 

J 

145 


BR=BR*l.E-05 

J 

146 


BTMETAsRTMETA*l .E-05 

J 

147 


0PH I aRpH I *1 . E"05 

J 

148 


B3=SURT(BR*BR*BTHETA»BTHETA*BPHI*BPHI) 

J 

149 


RETURN 

J 

150 


END 

J 

151 
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